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1.  INTRODUCTION 


The  transmission  of  standard  commercial  television  images  from  a  remotely  piloted  vehicle  (RPV) 
to  a  distant  observer  requires  a  large  bandwidth  channel  in  order  to  protect  the  video  information 
against  jamming  and  interference.  This  protection  can  be  gained  from  spread  spectrum  techniques.  In 
many  RPV  situations,  however,  the  available  channel  bandwidth  is  barely  sufficient  or  is  insufficient 
to  transmit  the  television  signal  directly.  In  these  cases  the  video  information  can  be  transmitted  only 
if  the  redundancy  present  in  the  original  images  is  greatly  reduced.  This  report  discusses  the  feasibility 
of  redundancy  removal  at  real  time  rates  with  small,  lightweight,  low  power  hardware  suitable  for  the 
ARPA  RPV  transmission.  It  is  desirable  to  be  able  to  vary  the  amount  of  redundancy  removed  to 
permit  the  observer  to  select  from  the  available  tradeoffs  of  resolution,  bandwidth,  scan  rate,  and 
jamming  environment. 

Image  redundancy  reduction  by  a  factor  of  ten  without  serious  image  degradation  has  been 
achieved  by  the  use  cf  linear  transformations  and  filtering  in  the  transform  domain.  Such  techniques 
are  the  only  methods  presently  known  to  obtain  such  large  redundancy  reduction  factors.  1  hey  are 
also  tne  only  methods  known  which  will  allow  the  redundancy  reduction  ratio  versus  degradation 
to  be  continuously  varied.  Kaar,  Hadamard,  Fourier,  and  Karhunen-Loeve  transforms  have  all  been 
successfully  used  in  this  way,  but  the  processing  has  been  accomplished  on  large  digital  computers  at 
much  slower  than  real  time  rates.  Examples  of  these  redundancy  reduced  photographs  are  given  in 
Appendix  A.  The  feasibility  that  the  required  transforms  can  be  performed  with  a  high  throughput 
processor  structure  which  utilizes  parallel  access  to  a  serially  shifting  data  stream  is  explored.  The 
processor  structure  may  be  implemented  by  snail,  lightweight,  low-power  hardware  using  transversal 
filters.  Other  methods  of  television  bandwidth  reduction  such  as  slow  frame  rate  and  image  sensor 
improvements  are  also  presented. 


1 


2.  VIDEO  DATA  SPECIFICATIONS 


PSYCHOPHYSICAL  STUDIES 

The  amount  of  information  that  an  operator  needs  to  effectively  operate  an  RPV  for  various 
missions  is  not  a  well-defined  quantity.  Rather  it  is  a  subject  of  some  controversy  in  the  literature^  *  1 
The  Aerospace  Medical  Research  Laboratory  (AMRL)  at  Wright  Patterson  /'ir  Force  Base  lias  been 
conducting  such  psychophysical  experiments  for  an  RPV  mission.  The  speed  and  altitude  of  the 
Wright  Patterson  RPV's  mission  differ  from  the  speed  and  altitude  of  the  ARPA  RPV.  The  Wright 
Patterson  studies  may  therefore  have  limited  applicability  to  the  ARPA  RPV  mission. 

However,  Wright  Patterson  has  sponsored  work  at  Virginia  Polytechnic  Institute  and  State  Uni¬ 
versity  in  order  to  determine  a  more  general  measure  of  video  image  quality.  This  work  has  met  with 
some  success  and  lias  been  reported  in  conference  proceedings^  ~  ^  and  in  a  final  report  to  the  Aero¬ 
space  Medic il  Research  Laboratory.  At  this  time  it  is  undergoing  review  by  AMRL  and  by  USC. 

In  the  absence  of  definitive  results  there  are  several  tentative  conclusions  that  have  been  reached. 
Since  the  ARPA  RPV  has  a  slow  speed,  low  altitude  mission,  and  since  the  human  operator  response 
time  is  approximately  one-quarter  of  a  second,  the  information  rate  (i.e. ,  frame  rate  for  television) 
need  only  be  approximately  four  frames  per  second.  Again,  since  the  ARPA  RPV  has  a  relatively  slow 
speed,  the  time  delay  of  the  information  for  control  and  reconnaissance  is  not  of  critical  importance. 
However,  it  is  important  in  the  target  designator  mission  and  should  always  be  less  than  the  human 
response  time,  preferably  much  less. 

The  image  quality  necessary  is  more  difficult  to  estimate.  The  resolution  of  home  quality  televi¬ 
sion  is  the  best  available  estimate  at  this  time.  The  Virginia  Polytechnic  stuuy  should  answer  this  in 
more  detail.  The  camera  itself  has  a  very  important  role  and  the  vidieon  is  not  the  best  camera  for 
RPV  missions.  The  Navy,  through  NAVLLLX,  is  developing  charge-coupled  device  cameras  which  will 
offer  a  better  alternative  to  the  vidieon. 

The  value  that  color  will  provide  to  an  RPV  operator  is  also  a  subject  of  controversy.  Studies 
have  been  made  with  contradictory  results.  This  subject  will  continue  to  be  investigated  The  Univer¬ 
sity  of  Southern  California  has  performed  independent  studies  in  this  area.  Their  work  will  be  detailed 
in  a  later  report.  If  it  is  determined  that  color  r,  necessary,  then  the  primary  implementation  diffi¬ 
culty  will  be  in  the  color  camera. 
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BANDWIDTH  COMPRESSION 

The  amount  of  bandwidth  required  to  digitally  transmit  television  images  is  determined  by  the 
bit  rate.  The  television  bit  rate  is  given  by 

B.R.  =  N*B*F  . 

where  N  is  the  number  of  picture  elements  (p  xels)  per  video  frame,  B  is  the  number  ot  quantization 
bits  per  picture  element,  and  F  is  the  frame  rate.  The  number  of  pixels  per  frame  is  a  measure  of 
resolution.  The  number  of  bits  per  pixel  is  a  measure  of  the  dynamic  range  of  the  sensor  and  the 
display.  The  frame  rate  for  broadcast  television  is  determined  by  the  avoidance  of  Bicker  in  the 
display.  Broadcast  television  lias  set  this  at  thirty  frames  per  second. 

Let  us  consider  a  channel  over  which  the  digitized  television  is  to  be  transmitted.  Let  the 
channel  bit  rate  be  a  given  quantity.  It  is,  in  most  applications,  a  number  fixed  by  considerations 
other  than  video  bandwidth.  It  is  determined  In  part  by  assumptions  of  white  noise  in  a  given  band¬ 
width  and  a  specified  bit  error  rate  which  implies  a  given  signal-to-noise  ratio  at  the  receiver. 

A  relative  compression  ratio  (in  dB)  can  then  be  defined  for  use  in  comparing  various  band¬ 
width  compression  schemes  as 

„„  /channel  bit  rate\ 

<  R  =  10  loe  [ — N*B?p — ) 

United  States  broadcast  television  has  thirty  frames  per  second,  525  lines  per  frame  with  480 
useable  lines  and  a  4:3  aspect  ratio.  Studio  quality  TV  is  stored  digitally  with  6  bits  per  picture 
element,  640  picture  elements  per  line,  and  480  lines  per  frame.  TL  data  rate  for  studio  quality 
television  is  55.3  megabits  per  second.  Let  us  postulate  a  20.0  megabit  channel  data  rate.  The 
relative  '..ompression  ratio  for  conventional  studio  quality  television  is  -4.4  dB;  that  is,  the  channel 
is  inadequate.  Home  quality  television  does  not  have  the  resolution  of  studio  quality.  The  Univer¬ 
sity  of  Southern  California  has  determined  that  a  sampled  image  of  256  by  256  with  6  bits  per  pixel 
is  comparable  to  home  quality  television.  This  has  a  data  rate  of  1 1 .8  megabits  per  second.  The 
C.R  for  home  quality  television  for  the  20.0  megabit  channel  is  2.3  db. 

Computer  studies  in  picture  transform  encoding  have  shown  that  the  transforms  allow  the 
average  number  of  bits  to  be  reduced  to  approximately  one  bit  per  picture  element  without  picture 
degradation.  Optimum  transform  encoding  has  reduced  tins  to  one-half  bit  per  pixel.  L  et  us  con¬ 
sider  the  C.R.  available  for  a  256  by  256  frame  at  4  frames  per  second  with  a  transform  encoder 
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capable  of  preservme  oicture  quality  at  one  bit  per  picture  element.  The  data  rate  is  2b2  kilobits 
per  second  and  the  C'.R.  for  a  20  n  egabit  channel  is  18.8  db. 

If  it  is  determined  that  the  resolution  of  home  television  is  adequate  for  RPV  missions,  then 
the  linear  transform  frame  encoding  plus  a  reduced  frame  rate  will  give  a  relative  compression  ratio 
of  18.8  db.  The  feasibility  of  a  real  time  hardware  implementation  of  this  bandwidth  reduction  will 
be  discussed  in  section  III. 

The  amount  of  AJ  available  on  the  redundancy  reduced  data  will  depend  on  the  subsequent 
coding  algorithms  used  to  encode  the  data.  A  coding  scheme  is  being  mves'igated  that  offers  an 
additional  10  db  of  AJ.  A  computer  emulation  of  the  system  has  been  programmed  which  will 
allow  the  testing  of  this  and  other  proposals.  This  emulation  is  described  in  Appendix  B. 

Signal  reconstruction  by  an  inverse  transform  at  the  ground  station  is  straightforward.  Addi¬ 
tional  processing  may  enhance  the  picture.  One  such  procedure  is  described  in  Appendix  C.  Recon¬ 
struction  enhancement  will  be  investigated  more  fully  in  the  next  phase  of  the  program.  Appendix 
A  contains  a  report  by  Habibi  and  Hershel  which  discusses  a  general  reconstruction  method. 


3.  TRANSFORM  ENCODING 


INTRODUCTION 

The  use  ot  unitary  transforms  for  image  encoding  has  been  evaluated  for  use  as  intraframe 
encoding  techniques^!  In  addition,  these  techniques  may  also  be  applied  to  interframe  and  multi- 
spectral  encoding.  However,  all  unitary  transformations  are  information  preserving  and  no  band¬ 
width  reduction  results  from  the  application  of  the  transform  to  the  image.  Instead,  the  transforms 
redistribute  the  variance  associated  with  each  picture  element  (pixel);  so  that  subsequent  to  the 
transform,  basis  restricting  operations  on  the  transform  coefficients  will  result  in  bandwidth  reduc¬ 
tion.  Upon  reconstruction  of  the  original  image  from  the  basis  restricted  transform  coeffieie  its,  a 
degraded  version  of  the  original  image  can  be  obtained.  Unfortunately,  the  interrelationship  be¬ 
tween  the  type  of  transform,  the  form  of  the  noninverf’ble  operation,  and  the  type  of  degra¬ 
dation  in  the  reconstructed  image  is  very  complicated  and  subjective.  The  universally  used  analytic 

criterion  of  the  mean-square-error  is  at  present,  the  best  technique  for  transform  comparison. 

For  the  particular  operation  of  basis  restriction  by  truncation,  a  particularly  simple 
interpretation  of  the  bandwidth  reduction  can  be  made.  The  transforms  may  be  viewed  as  a  vari¬ 
ance  redistributing  operation  that  approximately  decorrelates  the  transform  coefficients  while  trans¬ 
forming  the  variance  associated  with  each  picture  element  into  the  low-order  coefficients  of  the 
transform.  Under  the  assumption  that  each  set  of  picture  elements  can  be  considered  as  a  sample 
function  from  a  wide  sense  stationary  random  process  with  correlation  function  r^,  there  exists  an 
optimum  discrete  transformation,  the  Karhunen-Loeve  transformation,  which  totally  decorrelates 
the  transform  coefficients  and  maximally  compacts  the  variance  to  the  low-order  coefficients.  All 
other  transformations  can  be  compared  in  their  performance  by  comparing  their  transform  coeffi¬ 
cient  decorrelation  anu  viriance  compaction  with  this  optimum  transformation. 

This  intuitive  interpretation  can  be  made  rigorous  through  the  use  of  the  rate-distortion  criteri¬ 
on^!  (t  has  been  found  from  experience  that  the  closer  the  eigenvectors  of  the  transformation 
approximate  the  eigenvectors  of  the  optimum  Karhunen-Loeve  transformation,  the  greater  the  vari¬ 
ance  compaction  will  be  and  the  more  the  coefficients  can  be  truncated  while  maintaining  a  fixed 
rate  distortion  or  mean-square-error. 
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The  use  of  two-dimension  transforms  can  provide  improved  performance  over  the  use  of  trans¬ 
formations  on  a  line-by-line  basis^5 1  The  most  direct  approach  is  to  seek  a  two-dimensional  trans¬ 
form  which  simultaneously  decorrelates  the  transform  coefficients  and  compacts  the  variance  into  a 
corner  of  the  two-dimensional  transform  coefficient  space.  One  method  is  to  find  a  two- 
dimensional  transform  which  can  be  represented  as  the  product  of  a  transform  in  one  direction  and 
a  transform  in  the  other  direction.  Assuming  that  a  two-dimensional  picture  can  be  considered  as  a 
sample  function  from  a  random  process  with  two-dimensional  correlation  r(  r0  1  i.e.,  with  a 
correlation  coefficient  r(  in  direction  one  and  a  correlation  coefficient  r2  in  direction  two,  then  the 
optimum  discrete  transformation  is  the  successive  use  of  two  Karhunen-Loeve  transformations;  the 
first  with  parameter  r  j ,  and  the  second  with  parameter  r2. 

Another  variable  of  interest  in  transform  encoding  is  block  size.  For  a  one-dimensional  signal 
the  block  size  is  the  number  of  elements  of  the  transform  and  the  performance  of  the  transform 
improves  monotonieally  with  increasing  block  size.  For  two-dimensional  images,  transform  per¬ 
formance  also  increases  with  increased  number  of  elements  in  each  dimension  ot  the  transform. 
However,  two  dimensional  transforms  usually  require  intermediate  memory  to  store  the  transform 
coefficients  in  the  first  direction  while  the  transform  is  being  computed  in  the  second  direction. 

Also  of  interest  in  two-dimensional  transform  encoding  are  mixed  transforms,  e.g.,  one  hori¬ 
zontal  transform  and  one  vertical  transform.  Although  performance  increases  with  the  number  in 
each  direction  of  elements  in  the  transform,  performance  varies  with  the  particular  transform  chos¬ 
en.  However,  memory  requirements  tend  to  increase  linearly  with  the  number  of  elements  in  the 
second  transform  direction  since  all  of  the  coefficients  must  be  stored  from  the  first  transform.  The 
amount  of  intermediate  memory  may  be  minimized  by  the  use  of  a  small  block  size  for  the  image  in 
the  second  direction,  but  performance  may  not  be  optimized  by  this  choice.  The  choice  ot  a  mixed 
transform  thus  interacts  with  the  overall  system  design  and  the  available  memory  for  coefficient 
storage. 

KARHUNEN-LOEVE  TRANSFORMATION 

If  a  continuous  time  function  of  zero  mean  and  autocorrelation  function  rT=  e-0'7'  is  consid¬ 
ered  to  be  a  sample  function  from  a  wide-sense  stationary  random  process,  then  this  time  function 
can  be  expanded  by  the  Karhunen-Loeve  expansion^'  and  the  resulting  coefficients  will  be  uncor¬ 
related.  For  a  discrete  function  of  zero  mean  and  autocorrelation  function  RT  =  r  T',  which  may  be 
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considered  as  a  sample  function  from  a  first -order  Markov  process,  a  similar  discrete  Karhunen- 
Loeve  transformation  may  be  defined^  .  This  transformation  diagonalizes  the  covariance  matrix 
and  is  optimal  in  the  mean-square  error  sense  for  a  restricted  set  of  basic  functions  that  do  not  span 
the  complete  space. 

The  discrete  Karhunen-Loeve  expansion  is  given  by  [7]  for  the  case  N  =  2m  as 


2m  . 

G(k)  =  Y  — - — r  g(n)  sin  cc  [k  -(2m  +  l)/2]  +  rnijl 

2m  +  X2  1 

n- 1  n 


k  =  1,2,.  .  •  ,2m 


where 


X  = 


1  -H 


2 

-  2r  cos  wn  +  \ 


and  |oJn|  are  the  positive  roots  of 


tan  2mu>  = 


-(1  -r2)  sin  co 
(cos  cj  -  2r  +  r2  cos  w) 


Since  the  discrete  Karhunen-Lo£ve  expansion  involves  both  the  solution  of  a  transcendental 
equation  and  the  evaluation  of  the  autocorrelation  function  of  the  data  to  be  transformed,  real  time 
computation  of  this  transform  is  quite  complex.  However,  Habibi  and  Wintz  ^have  shown  that 
Karhunen-Loeve  transformations  using  approximate  autocorrelation  functions  are  satisfactory  for 
many  applications.  In  addition,  Pratt  at  USC  is  examining  the  use  of  corrected  Fourier  coefficients 
as  a  practical  way  of  computing  approximate  Karhunen-Loeve  expansions. 

The  preceding  considerations  suggest  that  the  Karhunen-Loeve  expansion  should  be  used  only 
for  small  block  sizes  and  that  under  these  conditions  either  precomputed  expansions  or  interpola¬ 
tion  may  be  used.  This  complexity  may  be  justified  for  luultispectral  encoding  or  for  reduction  of 
the  memory  required  for  a  mixed  transform. 


DISCRETE  FOURIER  TRANSFORM 

181 

Since  the  discrete  Fourier  transform  is  asymptotic  to  the  Karhunen-Loeve  transformation1  , 
even  though  the  basis  vectors  are  picture  independent,  the  Fourier  transform  represents  a  logical 
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choice  for  real  time  implementation.  The  Fourier  transform  exists  for  all  lengths  N,  The  basis 
vectors  are  complex  and  are  given  by 

4 £  =  e-j2rr  nk/N 

If  the  input  sequence  is  real,  then  the  Fourier  coefficients  will  be  conjugate  even.  If  the  input 
sequence  is  real  and  even  then  the  Fourier  coefficients  will  be  real  and  even.  If  the  input  sequence 
is  non-negative  and  symmetric,  then  the  Fourier  coefficients  will  be  an  autocorrelation  sequence. 

Many  methods  exist  for  the  computation  of  discrete  Fourier  coefficients.  The  Goertzel  algo¬ 
rithm  requires  a  number  of  computations  proportional  to  N2  but  can  be  used  for  all  lengths  N 
When  N  is  highly  composite  “fast”  transformations  can  be  used^l .  Thus,  if  N  is  of  the  form  2^, 
then  the  number  of  computations  can  be  made  proportional  to  Nq.  Although  “fast”  algorithms 
have  been  successfully  used  on  General  Purpose  Computers  they  are  too  slow  for  real  time  computa¬ 
tion  since  the  algorithm  iterates  q  times  before  achieving  a  solution.  This  problem  can  be  overcome 
by  the  use  of  q  processors  in  a  pipeline  architecture ^  *  ^1  ,  although  this  increases  the  complexity  of 
the  processor. 

A  linear  filter  implementation  also  exists  for  the  discrete  Fourier  transform  which  is  both 
easily  implemented  and  suitable  for  real  time  computation.  This  algorithm,  called  the  chirp-Z- 
transfornJ  *  *  1  is  based  on  the  substitution  nk  =  [n2  +  k2  -  (n  -  k)2] /2  and  can  be  used  for  any 
length  sequence  N.  The  transform  may  be  summarized  as  a  premultiplication  by  a  discrete  chiro, 
convolution  with  a  discrete  chirp  of  twice  the  length,  and  postmultiplication  by  a  discrete  chirp. 
This  transform  may  be  computed  with  either  acoustic  surface  wave  filters  or  charge  transfer  de- 
vices^ ' 2 1 . 

DISCRETE  COSINE  TRANSFORM 

Certain  properties  of  Foutier  transforms  should  be  observed  when  using  the  discrete  Fourier 
transform  (DFTT  The  DFT  is  the  Fourier  series  representation  of  periodically  extended  data  and  as 
such  has  convergence  properties  which  depend  on  the  input  data.  If  the  periodically  extend-'d  data 
is  discontinuous,  then  the  convergence  will  be  of  the  order  of  1/n;  if,  however,  the  periodically 
extended  data  is  continuous,  then  the  convergence  will  be  of  the  order  of  1  / n2.  1  be  cosine  trans¬ 
form  makes  the  data  appear  to  be  continuous. 
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A  discrete  cosine  transtorm  of  a  data  sequence  g(n),  n  =  0,  1 , .  .  .  ,  (m  -  1 )  can  be  def  ined  as 

N-l 

G(0)  =  :1/2/N  £  g(n) 
n=0 

N-l 

G(k)  =  2/N  J  g(n)cos[(2n+  1)  kir/2N]  k=l,2 . (N-l) 

n=0 

The  basis  vectors  are  a  class  of  discrete  Chebyshev  polynomials  which  are  real  and  are  given  by 

<j ^  =  cos  ( ( 2n  +  1 )  kjr/2Nl 

Ahmed  ^  ^  1  has  investigated  the  use  of  these  basis  vectors  as  substitutes  for  the  Karhunen- 
Loeve  basis  vectors  and  finds  that  they  are  superior  to  the  Fourier  basis  vectors  and  comparable  to 
the  Karhunen-Loeve  in  reducing  the  mean-square-error  in  basis  restricted  transformations  while 
maintaining  the  computation  simplicity  of  a  transformation  which  does  not  depend  on  the  picture 
statistics. 

A  somewhat  different  definition  of  the  cosine  transform  is  given  by  Cooley  et  alJ  along 
with  the  necessary  operations  required  to  compute  the  cosine  transform  with  the  Fast  Fourier 
Transform.  Flowever,  the  auxiliary  operations  are  somewhat  involved  and  the  length  of  the  trans¬ 
form  must  be  chosen  so  that  the  modified  data  is  correct  for  a  fast  Fourier  transform. 

In  order  to  take  advantage  of  subsequent  processing  algorithms,  it  may  be  desirable  to  have  the 
Fourier  coefficients  be  a  real  autocorrelation  sequence.  This  requires  that  the  data  sequence  be  ex¬ 
tended  so  that  it  is  symmetric  about  the  first  data  value  This  automatically  assures  that  any  subse¬ 
quent  periodic  extension  is  continuous  but  also  results  in  a  data  vector  which  has  an  odd  number  of 
values.  A  third  cosine  transform  can  be  defined  for  the  symmetrized  data  even  though  the  resulting 
sequence  length  is  odd.  This  transform  can  then  be  computed  with  the  chirp-Z  algorithm  or  direct¬ 
ly  by  the  modification  of  the  Fourier  transform: 


N-l 

G(k)  =  g(Q)  +  Re  ^  g(n)e-jkn^2N-1) 
n=  1 
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Even  though  these  symmetrized  transforms  conceptually  use  a  symmetrized  data  sequence, 
their  implementation  is  no  more  complicated  than  CZT  of  the  o.iginal  data  length  and  only  re¬ 
quires  a  change  in  the  CZT  reference  function.  The  performance  will  be  that  of  a  sequence  which 
has  been  extended  to  be  continuous  while  the  transform  will  be  real  and  half  of  u.i  autocorrelation 

function. 

The  third  cosine  transform  is  recommended  as  the  primary  transform  for  image  processing  in 
the  horizontal  direction  and  may  easily  be  used  on  the  video  in  real  time.  Its  use  in  the  vertical 
direction  must  be  compared  with  other  transforms,  however,  since  memory  must  be  employed  tor 
intermediate  storage  and  this  memory  increases  linearly  with  the  size  of  the  transform  ii.  the  vertical 

direction. 


SLANT  TRANSFORM 

In  ord«r  to  have  a  better  match  between  the  characteristics  of  television  images  and  the  basis 
vectors  of  the  orthogonal  transformations  used  to  transform  these  images,  Shibata  and 
Enumotol  1 5 1  introduced  a  transform  of  which  the  second  basis  vector  decreased  in  uniform  steps 
over  its  length.  Pratt,  Welch,  and  Chen1 161  have  developed  the  slant  transform  in  such  a  manner 
that  it  preserves  the  “sequent*  "  properties  of  the  Walsh-Hadamard  transform  as  well  as  maintains  a 

“fast”  algorithm  for  computation. 

Unfortunately,  the  transform  is  no  longer  binary  and  multipliers  must  be  used  in  the  computa¬ 
tion.  Thus,  the  transform  may  either  be  computed  digitally  or  by  means  of  a  transversal  filter  with 
as  many  weighting  networks  as  there  are  basis  vectors  in  the  transformation.  These  considerations 
appear  to  limit  the  applicability  of  the  slant  transform  to  applications  such  as  the  vertical  encoding 
of  television  images;  here  the  significant  improvement  in  performance  of  the  slant  transform  offsets 
lack,  at  present,  of  a  simple  real  time  computation  implementation. 

WALSH-HADAMARD  TRANSFORM 

The  Walsh-Hadamard  transform  in  one  dimension  has  basis  vectors  which  are  the  discrete  Walsh 
functions.  The  discrete  Walsh  functions  are  an  orthonormal  sequence  of  length  N  canonically  de¬ 
fined  for  N  =  2q,  q  a  positive  integer.  The  values  of  the  Walsh  functions  are  +1  and  -I .  The  “sequen¬ 
cy”  property  of  the  Walsh  functions  exists  when  the  basis  vectors  are  ordered  by  the  number  of  sign 


changes  which  occur  in  each  basis  vector:  the  number  of  sign  changes  increases  linearly  through  the 
non-negative  integers  trom  0  to  2^  -  1 .  This  corresponds  to  the  number  of  sign  changes  which  occur 
for  the  discre’O  Fourier  Transform  as  it  increases  in  frequency  front  DC  to  frequency  2^_I. 

The  canonical  Hadamard  matrix  of  length  8  is  of  the  form 

1111  1111  Sequency  0 


H(3)  = 


The  first  8  Walsh  functions  ordered  by  seq 


sequency  0 


aency  are  of  the  form 


One  factorization  of  H(3)  into  sparse  matrices  which  gives  a  “fast”  algorithm  similar  to  the 
fast  Fourier  transform l 1 7'  is  of  the  form 
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The  corresponding  signal  flowgraph  is  of  (he  form 
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which  is  observed  to  have  no  multiplication  other  than  +1 .  This  makes  possible  a  simple  hardware 
implementations  of  the  algorithm  and  its  computation  for  large  block  size. 

The  importance  of  large  block  size  becomes  significant  when  it  is  observed  that  the  one¬ 
dimensional  Walsh-Hadamard  transform  of  size  N2  is  the  complete  two-dimensional  Walsli- 
Hadamard  transform  of  size  N  if  the  input  vector  to  the  one-dimensional  transform  is  obtained  by 
appending  successive  rows  of  the  two  dimensional  data  together,  starting  with  the  first  row.  In 
terms  of  minimizing  the  auxiliary  memory,  this  tu’  -dimensional  property  may  be  significant.  How¬ 
ever  it  imposes  constraints  on  the  scanning  of  the  original  image  and  a  “pipeline”  implementation 
would  be  required  for  real  time  computation. 

Alternatively,  a  transversal  filter  implementation  of  the  Walsh-Hadamard  transform  is  possi- 
ble^  181 .  This  implementation  requires  that  the  number  of  weighting  networks  be  equal  to  the  size 
of  the  transform  which  becomes  large  if  a  one  dimensional  transform  is  used  for  two-dimensional 
processing. 

The  basis  vectors  of  length  8  in  “sequency”  order  and  shifted  for  real  time  implementation  by 
transversal  filtering  are  of  the  form 
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The  large  number  of  weighting  networks  can  be  eliminated  through  the  use  of  an  electronically 
variable  tap  weight  implementation.  If  the  basis  vectors  of  the  transformation  are  imagined  to  be  in 
“sequency”  order  and  successively  shifted  to  the  right  with  increasing  “sequency,”  then  a  tap  weight 
“sequency”  can  be  derived  that  will  successively  compute  each  Walsh-Hadam  -rd  transform  coeffi¬ 
cient  using  a  transversal  filter  with  2  N  -  1  taps.  The  time  behavior  of  a  15-tap  filter  is  shown  as  tap 
weights  for  an  electronically  variable  tap  weight  implementation  of  the  form 

_ _ _ _ _ ►  time 

Tap  1  10000000 

2  110  0  0  0  0  0 

3  1110  0  0  0  0 

4  11110  0  0  0 

5  1  1  T  1  1  0  0  0 

6  i  T  T  T  T  i  o  o 

7  i  T  T  T  T  T  l  o 

8  i  T  T  i  l  T  T  l 

9  o  T  l  l  i  i  i  T 

10  o  o  l  T  T  T  T  l 

n  o  o  o  T  T  l  T  T 

12  00001111 

13  o  o  o  o  oil  1 

14  00000011 

15  oooooool 
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HAAR  TRANSFORM 


The  H;iar  transform  in  one  dimension  has  basis  vectors  which  are  the  discrete  Haar  functions 
The  Haar  functions  are  an  or'honormal  sequence  {0”'  j  characterized  by  two  parameters  m  and  n 

1  or  a  sequence  of  length  N  =  2°  n  =  |  0,  1 . q  J  and  nt  =  { 1 .  2.  .  .  2n"  * }  .  A  transform  of 

length  N  =  2*  has  the  basis  vectors  in  the  form 


In  general  for  N  -  21' 
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I  The  Haar  transform  may  be  of  interest  in  image  encoding  since  it  is  a  generalized  differential 

encoding.  In  normal  differential  encoding,  N-l  first  differences  of  successive  data  points  ire  trans¬ 
mitted  as  well  as  the  first  data  point.  However,  if  an  error  occurs  in  any  of  the  differences  or  in  the 
first  data  point,  then  all  subsequent  reconstructed  values  arc  in  error  The  Haar  transform  tends  to 


overcome  tins  difficulty  by  transmitting  N/2  contiguous  differences  and  N/2  -  I  compound  differ¬ 
ences  as  well  as  the  sum  of  data  values. 

From  an  implementation  viewpoint,  the  Haar  transform  is  interesting  since  it  has  several  con¬ 
venient  hardware  implementations.  There  are  “fast,"  transversal,  and  recursive  filter  implementa¬ 
tions. 

The  "last"  implementation  of  the  Haar  transform  is  particularly  convenient  because  the 
original  matrix  has  many  zeroes.  That  is,  instead  of  requiring  a  number  of  operations  proportional 
to  N  log-)  N  the  number  of  operations  required  for  the  Haar  transform  is  proportional  to  N. 

Corresponding  simplifications  are  possible  in  a  transversal  filter  implementation  since  the 
number  of  weighting  networks  is  only  I  +  log2N  since  many  of  the  weighting  networks  may  be  used 
for  the  computation  of  more  than  one  Haar  coefficient 

The  use  of  the  Haar  transform  as  the  vertical  transfer  in  for  image  encoding  should  be  investi¬ 
gated  since  log-,  N  +  1  delay  line  memories  may  be  used  as  recirculating  integrators  for  the  recursive 
calculation  of  the  coefficients  of  the  transform.  This  transform  is  second  only  to  a  differential 
pulse  code  modulation  in  memory  requirement  as  a  vertical  transform  while  simultaneously  provid¬ 
ing  some  protection  against  error  propagation.  The  Haar  transform  does  not  have  a  "sequeney” 
property ,  although  it  may  be  considered  to  sample  the  input  waveform  at  progressively  coarser 
intervals^  *  C)1 . 
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4.  HARDWARE  IMPLEMENTATION  FEASIBILITY 


Laboratory  tests  ai  NUC  have  been  made  to  determine  the  feasibility  of  implementing  the 
vario  is  candidate  compression  systems.  A  computational  system  architecture  lor  performing  the 
two-dirnei  sional  image  transformation  which  meets  the  ARPA  requirements  ol  low  cost,  low  power, 
and  low  weight  has  been  determined .  This  architecture  is  described  in  more  detail  in  NUC  1 N  1 02b 
which  is  attached  as  Appendix  D. 

Two  areas  of  technology  have  been  investigated  in  detail:  surface  wave  devices  and  charge 
coupled  devices.  These  devices,  along  with  multipliers,  provide  a  means  of  implementing  a  real  time 
video  transform  processor.  The  block  size  which  can  be  built  for  a  two-dimensional  transform  is  a 
function  of  the  transform  technology  and  the  memory  technology.  The  memory  primarily  impacts 
the  choice  of  the  vertical  dimension  transform. 

A  surface  wave  device  presently  under  construction  at  NUC  has  been  designed  to  implement 
the  chirp-Z-algorithm  for  the  Fourier  transform.  The  overall  time  delay  of  the  filter  response  is 
related  to  television  scan  rates  and  has  been  chosen  as  1/4  of  53.5  ps.  It  was  judged  reasonable 
to  divide  this  into  64  sample  intervals,  which  led  to  the  choice  of  0.209  ps  as  the  sample  interval. 

Implementation  of  the  filter  at  bandpass  was  achieved  bv  further  subdivision  of  the  sample 
interval  into  an  integer  number  of  periods  that  correspond  to  a  frequency  in  the  range  25-30  MHz. 
The  chosen  integer  was  6,  and  the  corresponding  period  (frequency)  was  .0348  ps  (28.7  MHz). 

Each  sample  of  the  chirp  function  was  represented  by  a  group  of  three  finger-pairs  separated  by  a 
space  of  three  carrier  wavelengths. 

The  chirps  sweep  from  rero  to  about  2.4  MHz  to  zero  to  2.4  MHz  and  back  to  zero  again  (an 
inverted-W-FM).  They  are  represented  on  the  surface  wave  device  as  an  amplitude  modulated  series 
of  28  7  MHz  finger  groups,  each  group  being  weighted  by  the  amount  of  overlap  of  its  fingers.  One 
pair  of  sine  and  cosine  chirp  filters  is  interrogated  by  one  input  signal  and  the  other  pair  by  a  second 
input  signal.  The  four  outputs  are  combined  according  to  the  complex  CZT  algorithm. 

A  charge  transfei  device  (bucket  brigade)  transversal  filter  has  been  built  by  Texas  Instruments 
as  part  of  a  study  funded  by  Rome  Air  Development  Center.  It  is  a  200  stage  device  capable  of 
operating  at  around  0.2  MHz.  At  present  it  is  undergoing  evaluation  at  NUC.  It  appears  to  be 
capable  of  implementing  a  modiried  cosine  transform. 
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Th  *  extension  of  these  technologies  to  the  sizes  and  speeds  necessary  for  video  signal  process¬ 
ing  appears  to  be  straightforward.  The  preliminary  production  cost,  power,  and  weight  of  the  vari¬ 
ous  system  components  are  given  in  Table  1 . 

In  the  course  of  the  television  bandwidth  reduction  study,  the  television  sync  signal  was  identi¬ 
fied  as  a  major  area  which  had  not  been  thoroughly  considered.  It  is  desirable  to  provide  a  high  AJ 
margin  on  the  sync  signal  so  that  the  display  will  remain  stable  under  adverse  conditions.  Prelimi- 
naiy  feasibility  studies  have  indicated  that  the  RPV  status  could  be  included  on  the  sync  channel. 

Table  1.  Preliminary  System  Production  Costs,  Power,  Weight 


Development  Costs 

Production  Costs 

Transversal  filter 

5  1  2  Taps  CCD 

S50-75K 

$  5.00 

SWD 

in  house  (S50-75K) 

100.00 

9  months  -  12  months  development  time 

Power  <  1  watt 

Weight  <  1  lb 

Memory  for  vertical  transform 

16  X  256  CCD 

S50-75K 

S  15.00-20.00 

SWD 

S100K 

100.00 

9-1  2  mo.  development  time 

Power  <  2  watts 

Weight  <  1  lb 

TV  sync  +  data  encoder 

SWD 

S120K 

120.00 

Power  1 .2  watts 

Weight  <  1  lb 

TV  Transform  data  encoder 

LSI 

10K 

$  10.00-20.00 

Power  <  5  watts 

Weight  <  1  lb 
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5.  RECOMMENDATIONS  FOR  FOLLOW-ON  TASKS 


A  computational  system  architecture  has  been  determined  lur  performing  the  two-dimensional 
image  transformation  which  meets  the  ARPA  requirements  of  low  cost,  low  power,  and  low  weight 
The  recommended  transform  for  at  least  the  horizontal  scan  direction  is  the  cosine  transform  It 
has  performance  equivalent  to  the  optimal  transform.  The  choice  of  the  transfoim  in  the  verti¬ 
cal  direction  is  impacted  by  the  memory  available  in  the  RPV.  Continued  system  design  and  mem¬ 
ory  development  will  be  necessary  for  the  airborne  computational  system. 

The  following  tasks  are  recommended  for  follow-on  work- 

(a)  Continue  collaboration  with  the  University  of  Southern  California  to  develop  and  evalu¬ 
ate  various  transform  techniques. 

(b)  Continue  airborne  system  design. 

(c)  Determine  and  develop  necessary  memory  technology. 

(d)  Design  and  develop  transform  data  encoder  compatible  with  modem. 

(e)  Provide  television  sync  encoding  with  large  AJ  margin.  It  is  also  recommended  that  the 
down  link  data  be  included  in  the  television  sync  system. 

(0  Integrate  system  components  at  NUC. 

(g)  Simulate  the  channel  and  evaluate  the  AJ  of  the  video  system. 

(h)  Field  test  the  system  eighteen  months  after  receipt  of  ARPA  order. 

(i)  Monitor  continuing  work  at  Wright  Patterson  on  necessary  operator  requirements  for 
various  missions  with  Decision  Science. 

(j)  Design  and  develop  ground  station  display  and  postprocessing. 


18 


REPORT  SUMMARY 


This  report  details  the  results  of  the  first  phase  of  a  NUC  program  on  image  bandwidth  reduc¬ 
tion  for  application  to  the  ARl’A  RPV  problem  of  sending  television  images  ove,  spread  spectrum 
channels.  This  report  presents  primarily  a  study  and  feasibility  testing  phase. 

It  has  been  found  that  significant  bandwidth  reductions  are  available  by  means  of  linear  trans¬ 
formations  applied  to  the  image  and  filtering  in  the  transform  domain.  Such  techniques  are  the 
only  presently  known  methods  to  obtain  large  redundancy  reduction  ratios.  Bandwidth  reductions 
by  a  factor  of  6  are  easily  attainable  without  significant  picture  degradation.  Reductions  by  a  fac¬ 
tor  of  1  2  only  slightly  degrade  the  picture. 

The  image  redundancy  reduction  has  been  done  on  a  high  speed  general  purpose  computer  at 
the  University  of  Southern  California  Image  Processing  Institute.  The  Naval  Undersea  Center  has 
determined  that  the  linear  transform  techniques  developed  at  USC  and  NUC  can  be  implemented  in 
real  time  with  small  low  cost  hardwa'-'  suitable  for  the  ARPA  RPV. 

For  example,  the  Fourier  transform  may  be  implemented  in  this  manner  by  utilizing  the  chirp- 
Z-transform  algorithm  which  provides  an  exact  decomposition  of  a  discrete  Fourier  transform  into 
premultiplication  by  a  discrete  chirp,  correlation  with  a  discrete  chirp,  and  postmultiplication  by  a 
discrete  chirp.  The  hardware  required  to  implement  this  algorithm  is  thus  the  equivalent  of  a  com¬ 
plex  filter,  a  pair  of  complex  multipliers,  and  a  pair  of  complex  sequence  generators.  The  computa¬ 
tion  time  required  is  linear  in  the  desired  discrete  transform  length,  and  the  cost  is  nearly  linear  in 
the  length.  The  sequence  generators  may  be  read-only  memories,  and  the  filter  may  be  provided  by 
a  transversal  filter  device.  High  speed  acoustic  surface  wave  implementation  will  permit  a  two- 
dimensional  transform  to  be  realized  as  a  succession  of  one-dimensional  transforms.  When  high 
speed  CCD  devices  become  available  they  will  provide  direct  replacement  for  the  surface  wave 
filters. 

The  principal  result  of  this  preliminary  study  phase  of  the  program  is  the  recommendation  that 
a  hardware  development  and  implementation  phase  be  initiated.  A  flight  test  of  the  bandwidth 
reduction  system  is  feasible  eighteen  months  after  initiation  of  the  hardware  development  phase. 

An  overall  bandwidth  reduction  of  the  order  of  20  db  appears  to  be  feasible  without  significant 
picture  degradation. 
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APPENDIX  A 


use  PHOTOGRAPHS  AND  REPORTS 

This  appendix  contains  transformed  photographs  done  under  contract  N001  23-73-C-1507  at 
the  University  of  Southern  California.  The  second  part  of  Appendix  A  is  an  article  by  All  Habibi 
and  Ronald  S.  Hershel.  Ali  Habibi  is  tfe  principal  investigator  for  the  USC  contract.  A  more 
formal  USC  report  will  be  presented  by  the  Image  Processing  Institute  at  a  later  date. 
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Figure  9.  The  original  and  the  encoded  pictures  using  one¬ 
dimensional  transformations  and  DPCM  systems,  1  bit/pixel 
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Figure  10.  The  encoded  pictures  using  the  cascade  of  one 
dimensional  transformations  and  DPCM  systems,  2  bits/ 
pixel,  M  =  16. 


1 


1  Dim.  K.  L.  Transf.  8t  DPCM 
1  Bit/Pixel 


1  Dim.  Cosine  Transf.  &  DPCM 
1  Bit/Pixel 


Original 


1  Dim.  Cosine  Transf.  &  DPCM 
1  Bit/Pixel 


1  Dim.  K.  L.  Transf.  DPCM 
1  Bit/Pixel 


A  UNIFIED  REPRESENTATION  OF  DPCM  AND  TRANSFORM  CODING  SYSTEMS 

by 
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We  consider  a  transform  coding  system  that  uses  a  lower-triangular 
transformation  to  uncorrelate  the  data.  Based  on  this  transformation  we 
propose  a  generalized  DPCM  system  and  show  that  at  high  bit  rates  it  per¬ 
forms  almost  as  well  as  coding  by  the  method  of  principal  components 
(Karhumen-Loeve  transformation) . 

This  study  connects  the  transform  coding  system  to  the  DPCM  encoder 
by  showing  that  the  proposed  system  simplifies  to  a  standard  DPCM  encoder 
for  markov  data. 

I.  Introduction 

Two  important  classes  of  coding  systems  that  are  based  on  eliminating 
the  correlation  of  the  data  prior  to  quantization  and  subsequent  encoding 
are  the  differential  pulse-code  modulation  (DPCM)  and  the  transform  coding 
systems.  Both  classes  have  received  a  great  deal  of  attention  in  the 
recent  literature  and  have  been  used  successfully  in  coding  pictorial 
data. ^  ^  Historically  both  tec'wiques  were  developed  separately  and  are 

treated  individually  and  often  by  entirely  different  groups  of  researchers. 
The  only  established  link  between  the  two  systems  is  the  fact  that  both 
attempt  to  generate  a  set  of  uncorrelated  signals  prior  to  quantization  by 
a  memoryless  quantizer  or  a  set  of  quantizers.  The  Karhunen-Loeve  trans¬ 
formation,  also  known  as  the  method  of  principal  components  for  sampled 
data,  generates  the  uncorrelated  signal  by  using  an  orthogonal  transforma¬ 
tion  that  diagonalizes  the  autocovariance  matrix  of  the  data.  The  DPCM 
system,  on  the  other  hand,  is  designed  based  upon  modeling  the  data  by  a 
markov  process,  then  using  a  best  linear  prediction  to  obtain  the  set  of 
uncorrelated  signals. 


This  research  was  supported  by  the  Naval  Undersea  Center,  San  Diego  under  con¬ 
tract  N00123-73-C- 1507  and  by  the  Advanced  Research  Projects  Agency  of  the 
Department  of  Defense  and  was  monitored  by  the  Air  Force  Eastern  Test  Range 
under  contract  F08606-72-C-0008.  The  authors  are  with  the  Department  of 
Electrical  Engineering,  University  of  Southern  California,  Los  Angeles,  Ca.  90007 
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In  this  paper  we  will  consider  the  DPCM  system  and  show  that  it  is 
indeed  a  member  of  the  class  of  transform  coding  systems  where  a  lower 
triangular  matrix  is  used  to  transform  the  data  to  a  set  of  uncorrelated 
signals.  This  approach  will  provide  a  fresh  outlook  to  the  theory  of  the 
DPCM  system  realizing  that  the  DPCM  system  is  a  special  case  of  a  more 
general  coding  system  where  a  lower  triangular  matrix  is  used  to  transform 
the  data  to  a  set  of  uncorrelated  signals.  This  system,  unlike  coding  by 
unitary  transformations,  combines  the  operation  of  the  transformation  with 
that  of  the  quantizers.  The  entries  of  the  lower  triangular  operation  are 
related  to  the  covariance  matrix  of  the  data.  For  an  nth  order  markov 
process  the  triangular  operator  has  n  off-diagonal  bands  with  identical 
elements  within  each  band.  Thus,  for  this  special  case  the  transformation 
is  performed  recursively  and  is  identical  to  the  linear  predictor  in  an  nth 
order  DPCM  system. 


II.  Transformation  by  Lower-Triangular  Operators 

T 

Consider  an  N-dimensional  data  vector  X  “  (x, ,x_ , . . .  ,x.)  and  let 

12  b 

X  represent  a  sample  vector  from  an  ensemble  of  N-dimensional  zero  mean 
random  variables.  X  can  also  be  considered  as  a  vector  of  N  samples  {x^} 
that  results  from  sampling  a  continuous  random  process  uniformly  over  a 
finite  interval.  A  vector  Y  =  (y^ »y2» • • . »yn)  can  always  be  generated  from 
a  linear  combination  of  x . ' s  as 

l 


yl  -  X1 


1-1 


y.  «  x.  -  I  .x, 
3  j  k-i  10  k 


for  j  =  2, . . .N 


or  in  vector  form 


(1) 


Y  -  LX 


(2) 


where  L  is  a  unit  lower- triangular  matrix;  i.e., 
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L 
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9m 

9m 
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o  o 

o  o 

o  o 


(3) 


^N,N-1 1 


Denoting  the  covariance  matrices  of  X  and  Y  by  and  Cy  respectively,  (2) 
implies  that 


CY  *  L  CX  LT  W 

Cholesky  has  shown  that  for  every  symmetric  positive  definite  matrix  Cx 

there  exists  a  real  non-singular  lower-triangular  matrix  L  such  that  matrix 
T  (8) 

L  C  L  is  diagonal.  Martin  and  Wilkinson  have  considered  numerical 
algorithms  for  finding  L  and  Cy  and  have  developed  efficient  techniques 
requiring  only  N^/6  multiplications.  ^ 

The  fact  that  the  transformation  (2)  diagonalizes  the  covariance 
matrix  of  X  indicates  that  the  elements  of  Y  are  uncorrelated. 

Some  of  the  significant  properties  of  the  transformation  by  the  unit 
lower- triangular  operator  are: 

1)  The  unit  lower- triangular  operator  L  is  not  unitary;  thus  trans¬ 
formation  (2)  does  not  preserve  the  length  of  a  vector.  As  a  result  though 
the  determinant  of  the  covariance  matrix  of  X  is  invariant  under  this  trans¬ 
formation,  the  trace  of  the  covariance  matrix  is  not  invariant. 

2)  This  transformation  does  not  share  the  optimum  concentration  of 
energy  in  the  first  M  <_  N  components  of  Y  exhibited  by  the  method  of 
principal  components.  Indeed,  for  an  n^  order  markov  process  the  variances 
of  the  Y  components,  except  the  first  n  components,  are  all  equal. 


3)  Transformation  with  the  lower  triangular  operator  L  does  not 
require  a  transformation  delay. 
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Note  that  if  the  components  of  X  are  samples  from  an  n^  order 
markov  process  the  stochastic  linear  model  of  (1)  will  be 

n 

yj  ■  XJ  ■  Vj-k  J  -  . N  (5) 

where  x.  =  0  for  i  =  0,-1, -2,...  .  Then  the  operator  L  will  be  a  banded 

n 

matrix  of  n+1  bands;  i.e.. 


Transformation  of  an  N  vector  with  operator  L n  requires  less  than  nN 
multiplication  as  compared  to  N2/2  -  N  multiplications  needed  for  trans¬ 
formation  with  the  unit  lower-triangular  operator  L  in  its  general  form. 
Furthermore,  since  transforming  with  operator  requires  only  the  n  most 
recent  components  of  X,  this  transformation  can  be  performed  by  using  a 
feed-back  loop  identical  to  one  in  an  n^1  order  DPCM  system  to  perform  the 
transformation  recursively.  In  this  case  the  complexity  of  the  transforma¬ 
tion  is  independent  of  the  dimensionality  of  X  and  depends  only  on  the  order 
of  markov  process  n.  The  block  diagram  of  the  system  using  the  feed-back 

loop  for  operation  L  X  is  shown  in  Figure  1. 

n  3 


Coding  by  the  Method  of  Principal  Components 
Before  proceeding  with  the  lower  triangular  transformation  we  will 
briefly  review  coding  by  the  method  of  principal  components  and  block 
quantization.  The  block  diagram  of  this  system  is  shown  in  Figure  2.  Huang 
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and  Schultheiss '  considered  this  problem  first  and  proved  that  given  a 
total  number  of  binary  digits  N0  (a  finite  capacity  digital  channel)  the 
average  coding  error  is  minimum  if  vector  X  is  transformed  to  an  uncor¬ 
related  vector  Y  by  an  orthogonal  operator  A  and  Y  is  quantized  by  block 
quantization.  Matrix  A  is  an  operator  that  diagonalizes  the  covariance 
matrix  of  X  and  the  block  quantization  of  Y  involves  assigning  N0  binary 
digits  to  N  components  of  Y  proportional  to  the  logarithms  of  their 
variances  d^;  i.e., 

?  d. 

l  Jin  10  |c  1 1/N 

II  X  N 

where  |CX|  is  the  determinant  of  covariance  matrix  C  and  E  m  =  N0. 

X  i-1  1 

This  rule  was  obtained  by  observing  that  optimal  quantization  of  Gaussian  variables 

m . 

y,  using  a  quantizer  with  2  1  levels  introduces  a  quantization  error  q 
*  i 

such  that 

2  2  -m./2 

=  E{Vp  (10)  (8) 

Substituting  (7)  in  (8)  gives 

Etq-}  =  lCx!1/N  (1°)~9/2  =  A  (9) 

Equation  (9)  indicates  that  the  block  quantization  of  Y  results  in  an  equal 
quantization  error  A  for  all  components  of  Y.  Note  that  (9)  does  account 

for  the  inaccuracies  due  to  using  integers  for  nu  ,  thus  A  should  be  viewed 

as  the  lower  bound  of  the  quantization  error.  Since  each  quantizer  is 
independent  and  the  quant,  -ation  error  is  additive,  the  quantized  vector  Y* 


*  -2m, Jin  2 

Huang  and  Schultheiss  use  an  approximation  that  involves  2  1  ,  but 

further  studies  of  the  problem  have  shown  that  (8)  is  a  better 

approximation. 


ip 
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Y  =  AX  +  Q 


e(QQT}  =  Al 


where  I  is  the  identity  matrix.  This  gives  a  reconstructed 


vector  X  such 


X  +  A_1Q 


The  total  error  for  the  system  is  the  trace  of  the  covariance  matrix  of  the 
error  vector  X-X  ,  thus  from  (12)  the  average  coding  error  is 

eb  =  N  tr  =  A  (13) 


It  is  interesting  to  note  that  the  above  coding  system  not  only  minimizes 
the  coding  error  but  it  also  gives  an  uncorrelated  error  in  the  signal 

domain  which  is  a  desirable  property  in  many  applications  of  the  coding 
system. 

Note  that  if  all  components  of  Y  are  quantized  using  0  bits  for  each 
compote  it  the  total  coding  error  would  be  the  same  as  it  would  result  by 
coding  the  components  of  X  directly.  This  can  be  seen  by  noting  that  using  0 
bits  per  component  the  average  ceding  error  is 

‘s  t A''  *<»■'>  I*"1)1]  (14) 

where  the  components  of  the  quantization  error  Q  are  uncorrelated  with 
variances  equal  to  E{y?} (10) "6/2.  Therefore 

es  =  Tt  10"6/2  Cy(A_1)T]  =  io"0/2  o2  (15) 

where  a2  is  the  common  variance  of  {x.}  and  (15)  is  the  average  coding  error 

that  results  by  quantizing  the  components  of  X  vector  directly.  We  note  in 

passing  that  this  result  is  true  for  any  transformation  that  generates  an 

uncorrelated  Y  vector  including  the  transformation  by  the  lower-triangular 
operator  L. 
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IV.  A  Generalized  DPCM  System 

Transformation  of  vector  X  by  the  lower-triangular  operator  L  results 
in  a  vector  Y  where  the  components  of  Y  are  uncorrelated  and  in  general  have 
unequal  variances.  Components  of  Y  are  quantized  uring  the  block  quantiza¬ 
tion  technique  discussed  in  Section  III  and  are  transmitted.  Vector  X  can 
be  reconstructed  at  the  receiver,  within  some  level  of  degradation,  by 
operating  on  the  coded  vector  Y+Q  with  operator  L  \  The  average  coding 
2  . 

error  £  is 
Li 

c2-itr  [L-VV]  1161 

where  A  is  defined  by  (13).  The  readers  will  note  that  (16)  is  subject  to 
similar  inaccuracies  as  (18)  and  thus  it  also  should  be  viewed  as  a  lower 

bound  or  the  performance  of  this  coding  system. 

Comparison  of  (13)  and  (16)  show  that  using  lower-triangular  operator, 

a  non-unitary  matrix,  rather  than  the  method  of  principal  components  to  un¬ 
correlate  the  data  prior  to  its  block  quantization  results  in  an  inferior 
coding  system.  This  is  evident  from  the  fact  that  L  is  a  unit  lower- 
triangular  matrix,  thus 

trtlTVVlitr  I  <17» 

When  components  of  X  belong  to  a  first  order  markov  process  the 
block  diagram  of  the  above  coding  system  is  shown  in  Figure  3a.  Then  only 
two  quantizers  are  used  where  one  encodes  and  the  other  y2  through  YN<* 
since  y2  through  yN  have  identical  variances.  From  published  literature 
the  performance  of  this  encoder  improves  by  including  the  quantizer  in  the 
predictor  loop  as  shown  in  Figure  3b.  This  combines  the  operation  of  the 
quantizer  with  the  transformation  and  is  identical  to  a  DPCM  system  with 
the  stipulation  that  a  separate  quantizer  is  used  to  encode  the  first  com¬ 
ponent  of  the  differential  signal.  Naturally  the  effect  of  using  only  one 
quantizer  is  negligible  for  large  Nj  thus  two  systems  are  identical. 

Since  Lis  a  triangular  matrix,  the  operation  of  the  quantizers  can 
be  combined  with  the  transformation  to  give  a  generalized  DPCM  system.  We 
propose  a  combination  of  the  transformation  and  the  block  quantization, 
similar  to  er.e  in  the  DPCM  system,  and  show  that  the  performance  of  the 
generalized  DPCM  system  at  high  bit  rates  approaches  the  performance  of  the 
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encoder  using  the  method  of  principal  components  and  the  block  quantization. 
This  system  uses  N  /2  -  N  multiplications  to  transform  the  data  and  requires 
no  coding  delays.  It  simplifies  to  an  n^  order  DPCM  encoder  if  components 
of  data  X  belong  to  an  nth  order  markov  process. 

The  block  diagram  of  the  generalized  DPCM  system  is  shown  in  Figure 
4.  The  components  of  input  vector  x^  ,X2 '  *  * '  ,xn  are  °Pera*:e<^  uPcn  i-n  sequence 
at  times  t, ,t„,...,t  respectively.  Switches  S  ,  i  =  1,2,...N  are  closed  at 
t.  and  remain  open  at  all  other  times.  Variables  z^,  i  =  1,...,N  are 
generated  recursively  and  are  stored  in  separate  locations.  z^  through  z^ 
are  read  out  at  each  t,  and  are  multiplied  by  variables  through  C  ^  , 
respectively,  to  form  z  which  in  turn  is  stored  at  time  t^^.  Vector  C 


is  the  j 


.  th 


row  of  the  matrix  (L-I)  at  time  t,  for  j  =  1,2,...,N. 

3  * . 


Denoting  the  components  of  the  quantized  signals  by  (y^),  the  recur¬ 


sive  relation  for  Lz^}  from  the  block  diagram  is 


*  i"1 

z.  =  y.  +  £  H.  .z . 

1  1  j-i  1J  : 


for  i  =  2 , 3 , . . . ,N 


(18) 


and  zi  = 


furthermore 


and 


* 

Y  =  Y  +  Q 

Y  =Lx  +  (L-I)  Q 


(19) 

(20) 


where  Y  is  the  signal  at  the  input  of  the  quantizer,  Q  is  the  vector  of 

quantization  error.  These  equations  indicate  that  vector  Z  =  (z^,z^, . . . , 

z  )T  is  identical  to  the  reconstructed  vector  X  at  the  receiver.  Thus 
N 

from  (19) 


X  -  X*  =  Q  (2D 

and 

Y  =  Y  +  (L-I) Q  (22) 

where  Y  refers  to  the  uncorrelated  vector  LX.  From  (22)  and  realizing  the 
fact  that  E{y^q^}  =  0  for  j  =  l,2,...,i-l 
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i-1  l-l 

=  E{yJ}  +  l  Z  E{q.qj} 

j=l  k=l 


i  =  2 , . . .  ,N 


(23) 


and  E{y^}  =  E{y^}.  Note  that  {q± }  are  in  general  correlated  since  they 
result  from  quantizing  correlated  variates  {y^L  This  hampers  an  accurate 
analysis  of  this  system  however  assuming  fine  quantization  of  (y^)  the 
correlation  of  {q .  } are  negligible  and  if  block  quantization  of  ^y^)  Per 
formed  optimally,  using  the  variances  of  (it), .from  (8)  and  (23) 


E{qJ) 


dor""  |cxll/" 


1  -  (10) 


-0/2  /  ^Cx' 


1/N' 


*•  > 


.  ,N 


(24) 


i-1 

E  l 


ij 


and  E(q^)  =  A.  Equations  (23)  and  (24)  specify  variances  of  {yJ  in  terms 
of  the  covariance  matrix  of  the  data  and  the  bit  rate  0.  From  (21)  and 
(24)  the  average  coding  error  for  the  generalized  DPCM  system  is 


=  - 


N 

E 

i-2 


KiS 


A 

+  N 


(25) 


where  A,  defined  by  (9),  is  the  coding  error  by  the  met!  od  of  principal 
component.  Analogous  to  (9)  the  error  given  by  (25)  also  should  be  treated 
as  a  lower  bound  for  the  performance  of  the  generalized  DPCM  system. 

Figures  5  and  6  show  the  performance  of  this  system  as  given  by  (25)  and  the 
performance  of  the  optimum  system  (method  of  principal  components)  for  two 
non-markov  processes.  These  processes  are  defined  by  their  autocovariance 
functions 


*1 


(T)  = 


1  +  0.  2  T 


T  =  0,1,2,. 


(26) 


R2^  1  +  0.05  T 


T  —  0,1,2,. 


(27) 
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These  figures  also  show  the  coding  performances  of  PCM  and  the  encoder  that 
uses  transformation  by  operator  L  and  block  quantization  separately. 

For  the  special  case  where  {x^}  is  an  n  order  markov  process  the 
operator  L  is  replaced  by  L  of(l6).Then  the  elements  of  vector  C  remain 

toed  for  all  t„  except  for  ”  =  1,  2,  3 . a.  This  also  makes  the  variances 

of  the  components  of  the  transformed  vector,  again  except  the  first  n  com- 
ponents,  identical.  Thus  the  block  quantizer  will  consist  of  only  n+1 
quantizers.  This  reduces  the  generalized  encoder  to  an  n  order  DPCM 
system  with  the  stipulations  that;  first,  n+1  quantizers  are  used  to  encode 
the  differential  signal,  second  the  quantizers  are  designed  using  the 
variance  of  the  differential  signal  where  the  effect  of  the  quantization 
error  is  accounted  for.  Naturally  the  degradation  of  the  encoder  due  to 
using  only  one  quantizer  becomes  negligible  for  large  values  of  N. 

For  the  case  of  the  first  order  markov  process  the  coding  error  from 

(25)  for  large  values  of  N  is 


(io) “0/2  (1-aJ) 
I-  do) "0/2^ 


(28) 


The  numerator  of  (28)  is  the  coding  error  of  the  first  order  markov  process 
by  the  method  of  the  principal  components  for  large  values  of  N.  The  excess 
error  due  to  usr.g  an  improved  DPCM  rather  than  the  method  of  the  principal 
components  is  less  than  1.1%  for  =  0.9  and  0=2,  and  reduces  sharply  at 
higher  bit  rates. 


V.  Conclusions 

We  have  considered  a  transform  coding  system  that  combines  a  lower- 
triangular  transformation  with  a  block  quantizer  to  convert  a  set  of  sampled 
data  to  a  string  of  binary  digits  which  can  be  converted  back  into  a  replica 
of  the  data.  Separating  the  operation  of  the  quantizer  from  that  of  the 
transformation,  as  in  the  unitary  transform  encoders,  results  in  a  correlated 
coding  error  at  the  receiver  and  a  rather  inefficient  coding  system.  How¬ 
ever,  combining  the  operation  of  the  quantizers  with  the  transformation  of 
the  data  such  that  the  coding  error  at  the  receiver  is  uncorrelated  improves 
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the  coding  efficiency  of  the  system  significantly.  The  resultant  encoder 
is  only  slightly  suboptimum  to  the  system  that  uses  the  method  of  principal 
components  at  low  bit  rates.  At  high  bit  rates  the  performance  of  both 
systems  are  identical.  This  was  verified  for  a  number  of  processes  with 
different  autocovariance  functions. 

Transformation  by  the  lower-triangular  operator  involves  no  coding 
delay  and  requires  less  than  half  as  many  multiplications  as  the  method  of 
principal  components  requires. 

Combination  of  the  transformation  and  the  quantization  is  similar  to 
the  procedure  employed  in  the  DPCM  system.  Furthermore  the  proposed  encoder 
reduces  to  an  improved  form  of  DPCM  encoder  for  markov  processes,  thus  it  is 
referred  to  as  a  generalized  DPCM  encoder. 

The  fact  that  the  generalized  DPCM  system  reduces  to  a  DPCM  encoder 
for  markov  processes  establishes  a  link  between  the  DPCM  and  the  transform 
coding  systems  making  a  unifoed  approach  to  coding  by  these  two  techniques 
possible.  Furthermore  we  have  shown  that  the  theoretical  performance  of  a 
modified  DPCM  system  is  optimum  at  high  bit  rates.  At  low  bit  rates  it  is 
only  slightly  suboptimum  to  the  encoder  using  the  method  of  principal 
ponents. 
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This  report  describes  a  signal  processing  interpreter  which  was  developed 
as  part  of  the  AKPA  sponsored  program  for  image  transmission  via  spread  spectrum 
links.  The  signal  processing  interpreter  SPIN3  is  an  interactive  program  for 
use  at  a  time-sharing  demand  terminal.  It  provides  the  user  with  the  equivalent 
of  a  calculator  designed  to  perform  signal  processing  operations  and  provides 
the  software  equivalent  of  a  large  number  of  modules  for  breadboarding  a  complete 
signal  processing  system. 

This  interpreter  is  meant  to  be  used  by  engineers  and  scientists  who  are 
familiar  with  signal  processing,  but  who  may  have  no  knowledge  of  programming. 

The  signal  processing  interpreter  is  particularly  useful  for  the  rapid 
investigation  of  systems  whose  complexity  precludes  a  complete  analytic  study, 
and  whose  utilization  of  new  components  may  make  hardware  breadboarding 
undesirable  because  of  cost  and  procurement  time  limitations. 


INTRODUCTION 


The  signal  processing  interpreter  is  an  interactive  program  for  use  at  a 
time-sharing  demand  terminal.  It  provides  the  user  with  the  equivalent  of  a 
calculator  designed  to  perform  signal  processing  operations  -  the  software 
equivalent  of  a  large  collection  of  modules  for  breadboarding  a  corplete  signal 
processing  system.  Such  a  system  might  include  an  information  source,  a  source 
encoaer,  a  channel  encoder  or  modulator,  a  channel  with  noise,  jamming,  and 
multipath,  a  channel  decoder  or  demodulator,  and  an  information  decoder. 

Typically,  one  vrould  wish  to  compare  several  such  systems  operating  through 
each  of  several  channels. 

The  complexity  and  variety  of  the  systems  and  channels  will  frequently 
preclude  a  complete  analytic  comparison,  while  cost  and  procurement  time  limit 
the  flexibility  of  a  hardware  breadboard  of  a  system  utilizing  new  ccrtponents. 

The  t*io  critical  requirements  for  such  a  signal  processing  interpreter 
are:  (a)  It  must  include  the  most  frequently  needed  signal  processing  operations, 
with  a  provision  for  readily  adding  new  operations  as  the  need  arises;  and 
(b)  It  must  be  easy  to  use  by  engineers  and  scientists  who  are  familiar  with 
signal  processing,  but  who  may  have  no  knowledge  of  programming. 

The  signal  processing  interpreter  SPIN3  which  will  be  described  here  was 
developed  as  part  of  the  ARPA  sponsored  program  for  image  transmission  via 
spread  spectrum  links.  The  cotmand  set  is  meant  to  be  general  enough  for  a 
variety  of  signal  processing  problems.  In  addition,  special  conrands  are 
available  for  examining  and  processing  portions  of  pictures  in  order  to  facilitate 
the  comparison  of  proposed  image  transmission  systems. 
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HOW  TO  USE  THE  SIGNAL  PROCESSING  INTERPRETER 


The  signal  processing  interpreter  is  structured  around  an  operational  stack. 
Such  stacks  are  used  in  many  compilers,  in  the  Burroughs  B5000  computer,  and  in 
the  HP-35  calculator  [1,2,3],  Ihe  stack  structure  used  in  the  signal  processing 
interpreter  is  in  fact  identical  to  that  used  in  the  HP-35  except  for  the 
differences  necessitated  by  dealing  with  complex  vectors  rather  than  real 
scalars. 

The  stack  consists  principally  of  the  complex  vectors  X,Y,Z,  and  T.  There 
is  also  a  storage  vector  S  which  is  effected  only  by  the  cormand  STORE,  and 
which  effects  the  stack  only  through  the  cormand  RECALL.  Associated  with  each 
of  these  vector  registers  is  an  integer  NX,  NY,  NZ,  NT,  NS  which  indicates  the 
current  dimension  of  the  associated  vector. 

After  the  user  performs  the  start  up  procedure  described  in  Appendix  B, 
and  also  after  the  successful  completion  of  any  cormand,  the  interpreter  will 
solicit  a  command  by  displaying  the  message  "PLEASE  ENTER  NEXT  COMMAND" .  Ihe 
user  may  then  enter  any  cormand  (up  to  25  characters  enclosed  in  single  quotes) 
and  then  depress  the  carriage  control  key  to  send  the  cormand  to  the  computer. 

Che  of  the  five  following  outcomes  will  then  result: 

a)  The  cormand  will  be  recognized  and  performed,  and  a  new  cormand 
will  be  solicited. 

b)  The  cormand  will  be  recognized,  and  the  program  will  request 
the  user  to  enter  additional  data  from  the  keyboard. 

c)  The  cormand  will  be  recognized,  but  not  yet  implemented.  In 
this  case,  the  user  will  be  so  informed,  and  a  new  cormand  will  be  solicited. 

In  many  cases,  the  desired  cormand  may  be  replaced  by  a  short  sequence  of 
oarmands  which  have  already  been  implemented. 

d)  The  cormand  will  not  be  recognized;  the  user  will  be  so  informed; 
and  a  new  cormand  will  be  solicited. 

e)  The  attempted  cormand  is  not  in  the  form  of  an  expression  of  zero 
to  twenty-five  characters  enclosed  in  single  quotes.  This  will  produce  an  error 
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exit  fran  the  signal  processing  interpreter  and  an  Algol  error  massage  will  be 
displayed.  At  this  point  one  can  use  the  procedure  described  in  Appendix  B 
to  got  back  to  the  signal  processing  interpreter,  but  all  of  the  stack's  vector 
registers  will  be  reinitialized.  (It  is  planned  to  modify  the  program  in  the 
near  future  so  that  read  errors  will  only  result  in  an  attenpt  to  reread.) 

The  following  types  of  cortnands  are  available: 

a)  operations  on  a  single  vector 

b)  operations  on  a  pair  of  vectors 

c)  stack  manipulations 

>i)  exiting  from  the  signal  processing  interpreter. 

If  the  operation  has  a  vector  argument,  then  the  X  vector  is  used  as  the 
argument.  Similarly,  if  the  operation  has  a  single  vector  result,  the  result 
is  placed  in  X.  If  a  vector  function  of  a  vector  argument  is  evaluated,  the 
previous  value  of  the  X  vector  (i.c*.  the  argument)  is  destroyed  unless  seme 
special  action  is  taken  to  save  it,  such  as  storing  it  in  S.  Figure  1  shows 
the  stack  movements  for  this  case. 

old  new 

t  - »-T 

z  - — ►Z 

y  — - »~Y 

x  -  f  (x)  - *-X 

Figure  1.  Stack  movements  for  a  vector  function  of  a  vector  argument 

If  the  operation  is  an  output  operation,  such  as  printing  or  plotting,  the 
X  vector  is  used,  and  the  entire  stack  is  left  unchanged.  If  the  operation  is 
a  call  on  a  function  generator  such  as  TONE  or  CHIRP,  or  a  carmand  to  accept  a 
vector  input  from  the  keyboard,  the  vector  is  placed  in  the  X  vector  register, 
destroying  the  previous  contents  of  X. 

If  the  operation  evaluates  a  vector  function  of  a  pair  of  vectors,  then  X 
and  Y  are  used  as  the  arguments.  The  result  is  placed  in  X,  and  the  stack  is 
lowered.  The  previous  contents  of  X  and  Y  are  lost.  These  actions  are  shown 
in  Figure  2. 
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Figure  2.  Stack  movements  for  a  vector  function  of  too  vector  arguments 


For  most  functions  of  too  vector  arguments,  the  dimensions  of  the  too  argument 
vectors  are  first  automatically  made  equal  by  extending  the  shorter  one  with 
zeros  to  the  dimension  of  the  longer  one.  The  only  exception  at  present  is 
the  ootrtnand  APPEND  X  TO  Y,  which  leaves  a  resultant  of  dimension  NX  +  NY  in  X. 

«  * 

The  stack  operation  conmands  are  identical  to  those  of  the  HP-35,  except 
for  ROLL  UP,  which  is  an  exact  inverse  to  ROLL  DOWN.  Their  effects  are  shewn 
in  Figures  3-8. 

•  * 

0 

0 

Figure  3.  Stack  movements  for  the  cormand  'enter' 

0 

old  new 

t  - - *~T 


Y 
X 

Figure  4.  Stack  movements  for  the  cormand  'xy  interchange' 


old  new 
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PUW0  EXTENSIONS 


Several  extensions  of  the  signal  processing  interpreter  are  likely  in  the 
near  future.  In  order  of  probable  implementation  they  are: 

a)  Increasing  the  maximum  vector  dimension. 

b)  Implementing  the  carrmands  which  are  currently  recognized  but  not 
yet  implemented. 

c)  Providing  a  re-read  when  read  errors  occur. 

d)  Providing  additional  addressable  vectors  with  addressable  store 
and  reca  LI  ccnmands. 

e)  Providing  for  user-definable  compound  cormands  c.nd  synonyms  for 
present  ccrrmands. 

f)  Providing  for  ambiguity  function  calculation  and  the  use  of  more 
sophisticated  plot  urograms. 
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APPENDIX  A 
mm  LIST 


This  table  lists  all  commands  which  are  recognized  by  the  signal  processing 
interpreter,  even  if  the  ccrrmand  is  not  yet  implemented.  Contends  which  are 
recognized  but  not  yet  implemented  are  followed  by  (N)  . 

OPERATIONS  ON  A  SINGLE  VECTOR: 
fft 

reversal 
chirp  generator 
tone  generator 
czt 

conjugation 
input  vector 
print  vector 
plot  complex  vector 
picture  line 
set  x  dimension 
function  generator 
mirror 

sum  absolute  values 
pointwise  square  magnitude 
constant  function 
raise  to  the  power 
absolute  value 
plot (N) 

noise  generator (N) 
quantizer (N) 
histogram  (N) 

correlation  parameter  (N) 
hadamard  transform  (N) 


cosine  transform(N) 
kl  transform  (N) 
rectangular  to  polar (N) 
threshold  (N) 
binary  plot(N) 
shift  (N) 
print  max(N) 
scale  (N) 
autoscale  (N) 
exponential (N) 
modified  log(N) 

OPERATIONS  ON  A  PAIR  OF  VECTORS 1 
multiplication 
periodic  convolution 
subtract 
append  x  to  y 
periodic  correlation  (N) 
correlation (N) 
go] ay  sequence  generator  (N) 
add(N) 

STACK  OPERATIONS  : 
enter 
store 
recall 

xy  interchange 
clear (N) 
clear  x(N) 
roll  dcwn(N) 
roll  rq?(N) 


APPENDIX  B 

START  UP  AND  SIGN  OF  PROCEDURES 


This  section  assumes  that  the  signal  processing  interpreter  is  to  be  used 
from  a  demand  terminal  on  a  Univac  1108  operating  under  EXEC  VIII.  In  the 
examples  below,  user  entries  are  underlined  to  distinguish  them  from  the 
computer's  replies. 

A)  START  UP  PROCEDURE: 

UD108A  (enter  the  terminal's  SITEID) 

UNIVAC  1108  TIME/SHARING  EXEC  VERS  27.20.164  12  (the  system 
identifies  itself  to  the  user) 

@RUN  ABCD  (enter  a  short  run  card) 

@ASG,AZ  60*ARPA.  (assign  the  program  file) 

READY  (the  system  responds  when  the  file  is  available) 

@ADD,L  60*AKPA.BEGIN  (bring  in  the  remaining  needed  executive  control 
commands  automatically) 

@ASG,AZ  60*PICSET1.  (the  file  containing  the  USC  master  picture  set 
1  is  assigned) 

READY 

@USE  10.  , 60*PICSET1.  (the  name  10  is  associated  with  tte  picture 
set  file) 

READY 

@XQT  60*ARPA.SPIN3  (the  signal  processing  interpreter  is  called) 

B)  SIGN  OFF  PROCEDURE: 

'EXIT'  (the  user  exits  frcm  the  signal  processing  interpreter) 

EXIT  FROM  SIGNAL  PROCESSING  INTERPRETER 
EXIT  ALGOL-SIMULA  LIBRARY 

EXECUTION  TIME  2.2  SECONDS  -  LIBRARY  OF  SEPT  19,1969 
@FIN  (t)ie  user  terminates  the  run) 

(accounting  information  is  then  received  by  the  user) 


APPENDIX  C 

DETAILED  COTOND  DESCRIPTION 


This  appendix  only  describes  the  carman ds  which  are  currently  implemented. 
Cortmands  which  are  currently  recognized  but  not  yet  irplerrented  will  be 
documented  when  available. 

iih'  replaces  X  by  its  FFT  of  length  NX.  The  previous  contents  of  X  is 
destroyed.  The  current  X  dimension,  NX,  must  be  a  power  of  two. 

REVERSAL  X  is  replaced  by  X(0),  X(NX-l),  X(NX-2),  .  .  .  X'l),  i.e.  each  index 
is  replaced  by  its  negative  modulo  NX.  This  is  the  same  type  of  reversal  which 
is  produced  by  tv,j  forward  discrete  Fourier  transforms  in  succession. 

CHIRP  GENERATOR  generates  a  ccnplex  chirp  of  length  NX,  i.e.  X(k) 

for  k=0,  1,  .  .  .  NX-1.  The  previous  contents  of  the  X  vector  is  destroyed. 

TONE  GENERATOR  generates  uniformly  spaoed  sanples  of  a  sine  wave.  The  nimber 
of  sanples,  frequency,  and  sampling  increment  are  requested  of  the  user.  The 
resulting  sinusoid  is  placed  in  the  real  part  of  the  X  vector,  and  zero  is 
placed  in  the  imaginary  part  of  the  X  vector.  The  previous  contents  of  the  X 
vector  is  destroyed. 

CONJUGATION  replaces  the  X  vector  by  its  ccnplex  conjugate. 

INPUT  VECTOR  allows  the  user  to  manually  input  a  complex  vector.  First  he 
is  asked  for  the  dimension  of  the  vector  he  would  like  to  input.  He  then 
enters  each  term  in  the  form  REAL, IMAGINARY  followed  by  a  carriage  return . 

If  the  number  to  be  entered  is  purely  real,  the  imaginary  part  can  be  replaced 
by  an  extra  carriage  return.  The  previous  value  of  the  X  vector  is  destroyed, 
and  the  entered  vector  is  placed  in  X.  The  real  and  imaginary  parts  may  have  up 
to  5  digits  preceding  the  decimal  point  and  up  to  3  digits  following  the  decimal 
point. 

PRINT  VECTOR  prints  the  X  vector. 

PD3T  COMPLEX  VECTOR  normalizes  the  X  vector  to  a  maximum  absolute  value  of  1, 
prints  the  maximum  absolute  value  prior  to  no.malization,  prints  the  normalized 


X  vector,  and  simultaneously  plots  the  real  and  imaginary  parts  of  the  normalized 
X  vector. 

PICTURE  LINE  replaces  the  X  vector  by  a  section  of  a  line  of  picture  from  the 
USC  Image  Coding  Master  Set  No.l.  The  user  is  asked  to  specify  the  picture 
number  (1-9)  and  the  line  number  (1-256).  The  length  of  the  section  written 
into  X  is  the  smaller  or  256  or  NX. 

SET  X  Eu-MENSIQN  allows  the  user  to  reset  NX.  Immediately  following  the  oortmand, 
he  enters  an  integer  (1-512) .  If  the  integer  entered  is  smaller  than  the 
previous  value  of  NX,  it  replaces  the  previous  value  and  the  current  dimension 
of  tlie  X  vector  is  reduced.  If  the  integer  entered  is  larger  than  the  previous 
value  of  NX,  then  the  X  vector  is  filled  out  with  zeros  to  the  new  length,  and 
the  entered  integer  is  the  new  value  of  the  dimension,  NX, 

FUNCTION  GENERATOR  at  present  only  generates  a  linear  ranp,  i.e.  X(k)=k  for 
k=0,  1,  .  .  .  NX-1.  The  previous  contents  of  the  X  vector  is  destroyed. 

MIRROR  replaces  NX  by  2(NX)-1.  The  first  NX  entries  of  X  (i.e.  X(0),  .  .  . 
X(NX-l))  are  left  unchanged.  For  k=l,  .  .  .  NX-1  X(NX-l+k)  is  set  equal  to 
X  (NX-k) .  In  other  words ,  the  sequence  is  extended  to  length  2  (NX)  -1  so  as  to 
have  mirror  symmetry  (with  indices  modulo  2  (NX-1 ) )  about  the  zero  index  point . 

SUM  ABSOLUTE  VALUES  prints  out  the  sum  of  the  absolute  values  of  X(0) ,  .  .  . 
X(NX-l) . 

PQINTWISE  SQUARE  MAGNITUDE  replaces  each  of  X(0),  .  .  .  X(NX-l)  by  the  square 
of  its  absolute  value.  The  previous  contents  of  the  X  vector  is  destroyed. 

CONSTANT  FUNCTION  lets  the  user  enter  a  ccnplex  number  in  the  format 
REAL  IMAGINARY,  and  each  of  X(0),  .  .  .  X(NX-l)  is  set  equal  to  the  entered 
value.  Note  that  a  space  is  used  here  as  the  delimiter  between  tlie  real  and 
imaginary  parts.  The  previous  contents  of  the  X  vector  is  destroyed.  Other 
than  the  required  space,  the  input  is  free  format. 

RAISE  TO  THE  POWER  permits  the  user  to  enter  a  real  number  (with  decimal  point) , 
and  each  component  of  the  X  vector  is  raised  to  that  paver.  The  previous 
contents  of  the  X  vector  is  destroyed. 

ABSOLUTE  VALUE  replaces  each  element  of  the  X  vector  by  its  absolute  value, 
destroying  previous  contents  of  the  X  vector. 
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MULTIPLICATION  extends  the  shorter  of  the  X  and  Y  vectors  to  the  dimension 
of  the  longer  one.  Then  X  is  replaced  by  the  pointwise  product  of  the  two 
vectors,  and  the  stack  is  lowered.  The  original  values  of  the  X  and  Y  vectors 
are  destroyed. 

PERIODIC  OONVOLimCN  assumes  that  at  least  the  larger  of  NX, NY  is  a  power  of 
two.  The  shorter  vector  is  automatically  padded  out  with  zeros  to  the  larger 
dimension,  and  X  is  replaced  by  the  periodic  convolution.  The  stack  is  lowered, 
and  the  original  values  of  the  X  and  Y  vectors  are  destroyed. 

SUBTRACT  pads  out  the  shorter  of  the  X  and  Y  vectors  to  the  larger  of  the  two 
dimensions,  and  X  is  replaced  by  Y-X.  The  stack  is  lowered,  and  the  original 
values  of  the  X  and  Y  vectors  are  destroyed. 

APPEND  X  TO  Y  the  X  vector  is  replaced  by  X  appended  to  Y,  i.e.  by  Y (0) ,  .  .  . 

Y  (NY-1) ,  X(0) ,  .  .  .  X(NX-l) ,  and  the  new  value  of  NX  is  the  sum  of  the  previous 
X  and  Y  dimensions.  The  stack  is  lowered,  and  the  previous  values  of  the  X 
and  Y  vectors  are  lost. 


APPENDIX  D 
SOFTWARE  DESCRIPTION 


The  software  associated  with  the  signal  processing  interpreter  is  contained 
in  two  UNIVAC  1108  Fastrand  files;  a  program  file  called  60*ARPA. ,  and  a  data 
file  called  60*PICSET  1. 

The  data  file  contains  the  USC  Image  Coding  lister  Set  No.  1,  corrprising 
nine  pictures,  each  256  lines  long.  Each  line  consists  of  256  words,  and  each 
word  is  an  integer  in  the  range  0-255. 

The  program  file  contains  symbolic,  relocatable,  and  absolute  versions  of 
the  signal  processing  interpreter  SPIN3  and  its  associated  external  Algol 
procedures  and  Fortran  subroutines.  The  program  file  also  contains  assorted 

test  programs  and  file  manipulation  programs  which  were  used  in  the  course  of 
developing  the  software. 

The  Algol  programs  are  all  in  UNIVAC  1108  Algol,  and  the  Fortran  subroutines 
are  all  in  UNIVAC  1108  Fortran  V.  The  signal  processing  interpreter  and  the 
external  procedures  it  uses  are  shown  in  Table  1. 


PROGRAM 

NAME 

TYPE 

FUNCTION  PERFORMED 

SPIN3 

Algol 

signal  processing  interpreter,  main  program 

ATAB 

Algol 

procedure  for  camand  recognition  via  table  lookup 

FFTSP 

Fortran 

conplex  FFT  with  results  in  normal  order 

MULT 

Fortran 

subroutine  for  pointwise  multiplication  of  two  complex 
vectors 

REV 

Fortran 

subroutine  for  periodic  reversal  of  a  complex  vector 

CHIRP 

Fortran 

subroutine  to  generate  a  conplex  chirp 

TONE 

Fortran 

subroutine  to  generate  uniformly  spaced  sanples  of  a 
sinusoid 

STAR 

Fortran 

subroutine  to  replace  a  vector  by  its  conplex  conjugate 

PCONV 

Fortran 

subroutine  to  perform  the  periodic  convolution  of  two 
complex  vectors 

READV 

Fortran 

subroutine  to  let  user  manually  enter  a  conplex  vector 

D.l 


CPI0T 


Fortran 


Fortran 


Algol 


BEGIN 


subroutine  to  normalize,  print,  and  plot  a  complex 
vector 

subroutine  to  read  a  picture  line  from  a  suitably 
structured  mass  storage  file 

procedure  to  perform  a  discrete  Fourier  transform 
via  the  Chirp  Z-Transform  algorithm 

element  to  add  control  ccrrmands  to  the  run  stream  to 
get  ready  to  use  the  signal  processing  interpreter 
together  with  the  picture  data  file 


DESCRIPTION  OF  THE  fVVJOR  PROGRAMS : 

SPINS  tm.  u  . 

SPIN3  is  the  signal  processing  interpreter  main  program.  Upon  first  being 

entered,  the  program  prints  an  identification  message.  It  then  jumps  to  the 

statement  labeled  START .  The  successful  completion  of  an  action  by  the  main 

program  (possibly  including  a  call  to  an  external  procedure)  also  results  in 

a  jump  to  START.  _  , 

When  START  is  reached,  the  current  dimension  of  the  X  vector  is  prm  , 

together  with  a  request  to  the  user  to  enter  the  next  command.  At  this  point, 
control  is  transferred  to  the  external  Algol  procedure  ATAB,  which  reads  a 
string  variabie  (up  to  25  characters  enclosed  in  single  quotes)  and  returns  the 
variables  ICOM  and  CHECK  to  the  main  program.  The  boolean  variable  CHECK  is 
set  equal  to  true  if  the  string  is  raoognized  by  ATAB,  otherwise  it  is  set 
equal  to  false.  If  CHECK  is  false,  then  the  message  "COtlAND  NTT  EEOXHIZEn" 
is  printed,  and  the  program  again  jumps  to  START.  If  CHECK  is  true,  then  the 
integer  I(DM  is  used  to  determine  to  which  of  the  labels  LI,  L2,  etc.  control 
is  to  be  transferred.  If  the  oormand  is  recognized  by  ATAB,  but  not  yet 
implemented  in  SPM3,  a  jump  is  performed  to  the  statement  labelled  SOCN,  the 
user  is  notifies  that  the  comand  is  not  yet  implemented,-  and  the  program 

jumps  to  START. 

When  one  of  the  statements  labelled  LI,  L2,  etc.  is  reached,  the  requested 
action  is  performed  by  the  main  program  or  an  external  procedure  or  both.  In 
general,  stack  operations  and  short  amputations  are  performed  by  the  main 
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program,  while  relatively  conplicated  operations  are  performed  by  external 
procedures .  In  all  cases,  the  external  procedures  return  information  to  the 
nw-in  program  only  through  their  argument  lists.  While  this  method  pays  some 
penalty  in  address  computation  time,  it  greatly  facilitates  the  addition  of 
new  procedures,  changing  the  maximum  array  dimensions,  and  other  modifications 
as  additional  requirements  are  placed  on  the  signal  processing  interpreter. 

The  vectors,  X,  Y,  Z,  T,  S  currently  have  maximum  dimension  512  complex,  but 
this  parameter  could  be  easily  increased  if  desired. 

ATAB 

ATAB  is  an  Algol  procedure  which  reads  and  recognizes  signal  processing 
carman ds.  A  string  variable  S  is  read,  and  is  compared  with  each  of  the 
variables  in  the  string  array  T.  If  a  match  is  found,  the  integer  ICOM  is  set 
equal  to  the  corresponding  index  in  T,  and  the  boolean  variable  CHECK  is  set 
equal  to  true.  If  no  match  is  found,  CHECK  is  set  equal  to  false. 

FFTSP 

FETSP  is  a  Fortran  subroutine  to  perform  an  FET  on  a  complex  vector  having 
a  dimension  which  is  a  power  of  two.  It  is  a  minor  modification  of  Uhrich's 
FET  algorithm  [5],  The  main  changes  fran  Uhrich's  algorithm  are  that 
comnunication  is  only  through  the  argument  list,  and  that  N=2  is  specified 
in  the  call  rather  than  m. 

CZT 

CZT  is  a  straightforward  implementation  of  the  Chirp  Z-transform  algorithm 
described  by  Fabiner  et  al  [6] .  It  may  be  used  to  perform  the  discrete  Fourier 
transform  of  a  vector  of  arbitrary  dimension  (up  to  the  maximum  specified  in 
the  main  program  and  the  FET  subroutine) . 
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APPENDIX  E 
SAMPLE  RUN 


(USER  ENTRIES  ARE  UNDERLINED/  COWENTS  ARE  IN  PARENTHESES) 


UD108A  (user  enters  his  terminal's  siteid) 

UNIVAC  1108  TIME/SHARING  EXEC  VERS  27.20.225  01 

0RJJN  ABCD  (user  enters  a  short  run  card) 

DATE:  062273  TIME:  142232 

0ASG/AZ  60*ARPA.  (user  assigns  the  program  file) 

READY 

@ADD,L  60*ARPA.  BEGIN  (this  ccrrmand  adds  a  prestored  set  of  control  ccrmands 

to  the  run  stream  to  complete  the  initialization) 

@ASG,AZ  60*PICSETl.  (the  data  file  containing  USC  Master  Picture  Set  1  is 
READY  assigned) 

@USE  10,60*PICSET1.  (the  internal  filename  10.  is  identified  with  the  external 
READY  filename  60*PICSETl.) 

@XQT  60*ARPA.SPIN3  (the  signal  processing  interpreter  is  called) 

SIGNAL  PROCESSING  INTERPRETER,  VERSION  1 
NX= 

1 

PLEASE  ENTER  NEXT  COMMAND 

'FUNCTION  GENERATOR'  (requests  generation  of  a  linear  ramp) 

ENTER  DESIRED  FUNCTION  LENGTH 
8 

NX= 

8 

PIEASE  ENTER  NEXT  COMMAND 

'PRINT  VECTOR'  (prints  the  contents  of  the  X  vector) 

INDEX  COMPLEX  VALUE 


0 

.000 

.000 

1 

1.000 

.000 

2 

2.000 

.000 

3 

3.000 

.000 

4 

4.000 

.000 

5 

5.000 

.000 

6 

6.000 

.000 

7 

7.000 

.000 

FINISHED  PRINTING  VECTOR 

NX= 


PLEASE  ENTER  NEXT  COMMAND 

_|_STQRE'  (stores  the  contents  of  the  X  vector 


in  the  S  vector  register) 


PI  EASE  ENTER  NEXT  COMMAND 

^3?.'  (x  is  replaced  by  its  discrete  Fourier  transform) 
NX= 

8 

PLEASE  ENTER  NEXT  CQM-IAND 

'PRINT  VECTOR'  (prints  the  contents  of  the  X  vector) 


INDEX 


COMPLEX  VALUE 


0 

28.000 

.0C1 

1 

-4.000 

9.657 

2 

-4.000 

4.000 

3 

-4.000 

1.657 

4 

-4.000 

.000 

5 

-4.000 

-1.657 

6 

-4.000 

-4.000 

7 

-4.000 

-9.657 

FINISHED  PRINTING  VECTOR 

NX= 


8 

PIEASE  ENTER  NEXT  COMMAND 

(recans  S  into  X,  so  X  new  contains  the  linear  ramp) 


8 


PLEASE  ENTER  NEXT  COMMAND 


'MIRROR' 

NX= 


(generates  a  symmetric  function  of  length  2(NX)-1  by  appending 
to  X  its  mirror  image  about  the  zero  index  point) 


15 


PLEASE  ENTER  NEXT  COMMAND 

'PRINT  VECTOR'  (prints  the  contents  of  the  X  vector) 
INDEX  COMPLEX  VALUE 


0 

.000 

.000 

1 

1.000 

.000 

2 

2.000 

.000 

3 

3.000 

.000 

4 

4.000 

.000 

5 

5.000 

.000 

6 

6.000 

.000 

7 

7.000 

.000 

8 

7.000 

.000 

9 

6.000 

.000 

10 

5.000 

.000 

11 

4.000 

.000 

12 

3.000 

.000 

13 

2.000 

.000 

14 

1.000 

.000 

FINISHED  PRINTING  VECTOR 

NX= 


E.2 


I 

I 


15 

PLEASE  ENTER  NEXT  COMMAND 

'CZT'  (replace  X  by  its  discrete  Fourier  transform,  corrputed  via  the 
NX=  Chirp  Z -Trans  form  Algorithm) 


15 

PLEASE  ENTER  NEXT  COMMAND 

'PRINT  VECTOR1  (print  the  contents  of  the  X  vector) 
INDEX  COMPLEX  VALUE 


0 

56.000 

.000 

1 

-22.881 

.000 

2 

-.261 

.000 

3 

-2.618 

.000 

4 

-.300 

.000 

5 

-1.000 

.000 

6 

-.382 

.000 

7 

-.558 

.000 

8 

-.558 

.000 

9 

-.382 

.000 

10 

-1.000 

.000 

11 

-.300 

.000 

12 

-2.618 

.000 

13 

-.261 

.000 

14 

-22.880 

.001 

FINISHED  PRINTING  VECTOR 
NX= 


15 

PLEASE  ENTER  NEXT  COMMAND 


'PICT  COMPLEX  VECTOR1  (normalize  X  to  a  maximum  absolute  value  of  1.0,  print  and  plot) 
REAL  PLOT  USES  R,  IMAGINARY  PLOT  USES  I  AS  CHARACTER 


i* 

MAX  ABS  VALUE 

IN  ARRAY= 

56.00000 

1 

1.00000 

.00000 

2 

-.40859 

.00000 

mt 

3 

-.00467 

.00000 

4 

-.04675 

.00000 

5 

-.00535 

.00000 

6 

-.01786 

.00000 

m 

7 

-.00682 

.00000 

m  m 

8 

-.00997 

.00000 

1 

9 

-. 00997 

.00000 

*W 

10 

-.00682 

.00000 

11 

-.01786 

.00000 

12 

-.00535 

.00000 

ak* 

13 

-.04675 

.00000 

14 

-.00467 

.00000 

T 

15 

-.40858 

.00002 

1 
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NX= 


I 

P.  I 

I 
RE 
I 
PI 
I 
I 
I 
I 
RI 
I 
RI 
I 

R  I 


R 


15 

PLEASE  ENTER  NEXT  COMMAND 

'EXIT'  (exit  frcrn  the  signal  processing  interpreter) 

EXIT  FROM  SIGNAL  PROCESSING  INTERPRETER 
EXIT  ALGOL-SIMLJLA  LIBRARY 

EXECUTION  TIME  000000.918  SECONDS  -  LIERARY  OF  SEPT  19,  1969 
@FIN  (terminate  the  run  and  call  the  accounting  routines) 


! 


! 


1 


i 


RLJNID: 

ABCD 

ACCOUNT:  XXXXXXXX 

PROJECT:  XXXXXXXXXXXX 

*1 

CPU  TIME 

00:00:0; 

.102  @  $.  300/SEC 

= 

$ 

.33 

.1 

CORE  SECS 

81 

@  $. 300/CBS 

= 

$ 

.24 

] 

I/O  REFS 

51 

@  $.  300/CRF 

= 

$ 

.15 

WDS  XFD(E2) 

902 

@  $. 000/WRD 

= 

$ 

.00 

CARDS  IN 

21 

@  $. 001/CRD 

= 

$ 

.02 

j 

PAGES 

5 

@  $.040/PGE 

= 

$ 

.20 

! 

.  i 

TOTAL  FUN  GOST  (TENTATIVE) 

= 

$ 

.94 

1 

INITIATION  TIME:  14:22:31-JUN  22,1973  VERSION:  27.20.225  01 


****T,TMR  INACTIVE**** 
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APPENDIX  C 


APPLICATION  OF  MAXIMUM  ENTROPY  ESTIMATION  TO  IMAGE  TRANSMISSION 

The  intensity  of  an  image  line  may  be  represented  by  a  nonnegative  function  Hx),  on  the 
interval  [0,  LI  The  Fourier  transform  of  I  is  nonnegative  definite,  and  hence  may  be  regarded  as 
an  autocorrelation  function.  Reconstruction  ot  the  image  line  from  a  truncated  version  ot  the 
Fourier  transform  of  the  image  is  thus  mathematically  equivalent  to  the  determination  of  a  spectral 
density  function  from  a  truncated  autocorrelation  function.  If  desired,  l(x)  may  be  replaced  by  its 
mirror  symmetrized  version  g(x)  =  l(|x|).  The  Fourier  transform  of  g  is  a  real  autocorrelation 
function. 

For  wide-sense  stationary  random  processes  with  all-pole  spectra,  the  method  maximum  en¬ 
tropy  provides  an  accurate  means  to  extrapolate  a  truncated  autocorrelation  function  and  to  use  the 
extrapolated  autocorrelation  function  to  torm  a  spectral  estimator  having  much  finer  resolution 
than  that  of  conventional  (Blackman-Tukey  type)  estimators  that  use  the  truncated  autocorrela¬ 
tion  function. 

We  propose  to  address  three  problems  in  order  to  determine  the  applicability  ot  maximum 
entropy  estimation  to  reduced  redundancy  image  transmission: 

a.  Applicability  of  the  all-pole  spectrum  model. 

b.  Feasibility  of  a  real-time  computational  algorithm. 

c.  Accuracy  requirements  for  hardware  to  implement  a  real  time  computational  algorithm. 
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INTRODUCTION 

The  time  requ ired  to  perform  signal  processing  computations  is  determined  largely 
by  the  time  required  to  perform  linear  transformations  on  data  vectors.  The  choice  of 
algorithm  and  implementation  structure  strongly  attects  both  the  processing  time  and  the 
implementation  complexity  via  the  required  number  ot  multiplication  times  and  via  the  data 
storage  and  transfer  structure, 

Any  linear  transformation  that  takes  a  sequence  of  N  data  points  into  a  sequence  ot 
N  transform  points  may  be  regarded  as  a  multiplication  ot  a  vector  by  an  N  X  N  matrix  A 
direct  implementation  that  uses  a  single  multiplier  requires  N“  multiplication  times  and  N~ 
words  of  storage.  The  Fast  Fourier  Transform  (FFT)  requires  a  number  of  multiplications 
proportional  to  N  logT  N. 

This  paper  discusses  a  numbe,  of  transform  implementations  that  have  a  simple 
serial  access  data  flow  and  a  computation  time  proportional  to  N.  These  transform  imple¬ 
mentations  which  consist  ot  two  complex  multipliers  and  a  complex  convolver  as  a  module, 
provide  an  effective  method  by  which  the  equivalent  of  a  large  number  of  multipliers  can  be 
used  in  a  processing  architecture  with  a  simple  data  How.  The  transform  structure  discussed 
here  may  be  used  for  different  types  of  transforms  including  the  Discrete  Fourier  Trans¬ 
form  (DFT)  by  a  change  of  the  constants  provided  in  a  local  store.  If  a  crossconvolver  is  used 
in  this  structure,  it  may  be  used  also  to  perform  direct  high  .peed  linear  filtering  and 
crosscorrelation. 

STRUCTURE  OF  TRANSFORMATIONS 

Several  types  of  elementary  linear  transformations  on  data  vectors  may  be  viewed  as 
building  blocks  for  the  implementation  of  more  general  linear  translormations.  Among 
these  elementary  transformations  are  multiplication  by  a  diagonal  matrix  (Fig.  1 ),  convolu¬ 
tion  (Fig.  2),  and  circular  convolution  (Fig.  3).  A  transformation  of  particular  importance 
for  signal  processing  is  the  discrete  Fourier  transform.  The  chirp  Z-transform  (CZTt 
algorithm  [  1 1  decomposes  the  DFT  into  a  premultiplication  by  a  discrete  chirp,  a  convolu¬ 
tion  with  a  discrete  chirp,  and  a  postmultiplication  by  a  discrete  chirp  (Fig.  4),  since 

N- 1  it  o  \ 

G  ,  V  e-i27mmgm  x  e-irrn2/N  >  e*n-m^/N  (e-nrm-/N  gJ 
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Fig.  1 .  Multiplication  by  a  diagonal  matrix 


If  the  reference  functions  used  in  the  convolver  and  the  multipliers  are  varied,  the  same 
architecture  used  to  implement  the  CZT  may  be  used  to  implement  a  wider  class  of  linear 
transforms  (Fig.  5). 

A  linear  transform  may  be  implemented  by  this  architecture  if  and  only  if  its  kernel 
decomposes  appropriately.  For  a  continuous  transform 


OO 

G(u)  =  !  k(u,  t)g(t)  dt 

-  OO 

the  required  decomposition  is 


k(u,  t)  ~  a(t)  b(u-t)  c(u) 


if  a  convolution  is  to  be  used  as  the  intermediate  operation.  Similarly,  if  a  correlation 
is  to  be  used  as  the  intermediate  operation,  the  required  decomposition  is 
k(u,  t)  =  a(t)  b(t  +  u)  c(u ). 

If  the  kernel  is  smooth  and  nonvanishing,  a  necessary  and  sufficient  condition  for 
decomposition  is 


d  log  k  5  log  k  _  £  /  dk  d_k  \ 
dt  3u  “k'dtou/ 


Any  function  of  t  plus  any  function  of  u 
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An  alternative  necessary  and  sufficient  condition  is 
33  log  k  +  33  log  k  _  Q 

a  t“  au  at  3u3 

where  the  plus  sign  corresponds  to  the  convolution  decomposition  and  the  minus  sign  corre¬ 
sponds  to  the  correlation  decomposition. 

For  a  discrete  transformation  of  the  form 

^  ^m,n  8n 

the  corresponding  required  factorization  'or  a  convolution  decomposition  is 
^m.n  “  an*3m-ncm 

where  the  necessary  and  sufficient  condition  for  factorization  (assuming  that  the  kernel  is 
nonvanishing)  is 

km+i  n  -  km  n-  1  =  any  function  of  m  pit’s  any  function  of  n 

where 

^m,n  ~  'og  ^mm 

An  alternative  necessary  and  sufficient  condi. ion  is 
wm,n  =  0 

where 

km,n  “  *og  km,n 
um,n  =  km+l,n  -  km,n-l 

vm,n  =  um+l,n  ~  um,n 
wm,n  =  vm,n+l  ~  vm,n 

If  the  discrete  transform  is  to  be  applied  to  a  data  vector  of  length  N,  then  a  1  sub¬ 
scripts  are  to  be  interpreted  modulo  N. 

A  listing  of  some  of  the  more  useful  decomposable  transforms  is  shown  on  the 
following  page. 
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Decomposable  Transforms 


TRANSFORM 

KERNEL 

FOURIER 

e-itu 

LAPLACE 

e_tu 

HILBERT 

_1_ 

t-u 

FOURIER  COSINE* 

cos  tu 

FOURIER  SINE* 

sin  tu 

•NOT  STRICTLY  DECOMPOSABLE  BUT 
REALIZABLE  WITH  SAME  HARDWARE 
AS  DECOMPOSABLE  TRANSFORMS 

COMPLEX  ARITHMETIC 

Complex  arithmetic  operations  applied  to  sampled  data  may  be  implemented  >n 
parallel  or  interleaved  form.  The  parallel  implementations  (Figs.  6  and  7)  follow  directly 

aR 


9  •  a  =  (gR  +  igt)  (aR  +  ia,)  =  (gR  aR  -  g,  a,)  +  i  (gR  a,  +  g,  aR) 

Fig.  6.  Complex  Arithmetic  Multiplier  (parallel  implementation) 
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SO 


b 


v 


IMAGINARY  OUTPUT 


g*b  ■  (gR  +  igj)  ★  (bR  +  ib|)  =  (gR ★  bR  -  g(  ★  b|)  +  i  (gR  *b,  +  g,  *bR) 

Fig.  7.  Complex  Arithmetic  Convolver  (parallel  implementation) 

from  the  definition  of  a  complex  multiplication  and  a  complex  convolution  The  operation 
of  the  complex  convolver  with  interleaved  arithmetic  (Fig.  8)  may  be  understood  by  an 
examination  of  the  operations  on  the  data  samples  in  the  even  and  odd  positions  in  the  data 
vector. 

We  wish  to  evaluate 

g*  b  =  (gR  +  i  g|)  *  (bR  +  i  bj)  =  (gR  *  bR  -  g(  *  bj)  +  i  (gR  *  b(  +  g(  ★  bR) 


Fig.  8.  Complex  Convolver  with  Interleaved  Arithmetic 


so  the  desired  real  term  is 


^  bn,Rgm-n,R  ^bn,lgm-n,l 


and  the  desired  imaginary  term  is 


Let 


Z  bn 

,1  gm-n,R 

n 

x2k  = 

gk.R 

x2k+l 

gk,l 

y2k  = 

bk,I 

y2k+l 

=  bk,R 

At  an  even  shift  position,  say  m  =  2p, 

X  xk  V 2p-  k  =  J  X2s  y2p- 2s  =  X  gs.R  bP- s,I 
k  even  s  s 

^  xk  y2p-k  "  ^  x2s+l  y2p-2s-  1  “  ^  g;  1  bp-s-  1  ,R 
k  odd  s  s 

and  the  two  terms  must  be  combined  with  a  relative  delay.  At  an  odd  shift  position,  say 
m  =  2p+l 


,R 


X  xky2p+l-k  J  x2sy2p+l-2s  "  ^  8S,R  bp-s,I 
k  even  s  s 

X  xky2p+l-k  =  ^  x2s+l  y2p+l-(2s+l)  =  ^  X2s+1  y2(p-s) 
k  odd  s  s 

~  X  gs,l  bp-S,l 
s 

and  the  two  terms  must  be  combined  with  zero  relative  delay. 


A  DISCRETE  FOURIER  TRANSFORM  IMPLEMENTATION 

Let  us  examine  the  discrete  Fourier  transform  decomposed  via  the  CZT  algorithm, 
and  implemented  with  the  parallel  form  of  the  complex  arithmetic. 
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N- 


N-l 


q  =  V  e“i27rmn/Ng  _  e  i7rm“/N  ^  e+i)r(m- n)“/N  ^e-  bm^/N  g  j 
n=0  n=0 

The  complex  operations  that  we  wish  to  implement  are  shown  in  Fig.  5  where 
i7rn-/N 


an  cn  e 


tor  n  =  0.  .  .  .  N  1 


bm  =  e+i™'/N 
m 


for 


r  m  =  0,  .  .  . 

Ini  =  0,  .  .  . 


ni  =  0,  .  .  .  2N-  2  if  correlation  is  used 

N-  1  if  circular  correlation  is  used 


When  the  parallel  implementation  of  the  complex  arithmetic  is  used,  the  DFT  archi¬ 
tecture  of  Fig.  9  results.  If  the  data  vector  dimension  N  is  even,  then  bm+;q  =  bm,  and  a 
circular  convolution  of  length  N  may  be  used  in  place  of  the  convolution  of  length  2N-  1 . 


CONVOLUTION  IMPLEMENTATION  WITH  TRANSVERSAL  FILTERS 

The  convolution  operation  may  be  implemented  in  highly  parallel  form  as  a  trans¬ 
versal  filter.  Such  a  filter  is  shown  in  Fig.  3  with  the  output  of  the  filter  connected  buck  to 
the  input  so  that  it  may  be  used  for  circular  convolution.  Various  acoustic  transversal  filters 
(21  have  been  described  that  are  suitable  for  use  in  a  fixed  transform  implementation  if  the 
signal  processing  can  proceed  without  interruption  during  each  block  of  the  input  signal.  If 
interruption  of  the  signal  processing  is  required,  nonacoustic  transversal  filters  may  be  em¬ 
ployed  [3] .  The  three  types  of  filters  which  will  be  examined  here  employ,  respectively, 
acoustic  surface  waves,  cha  ge  transfer  phenomena,  and  signal  propagation  in  digital  shift 
registers. 

Acoustic  Surface  Waves 

For  signal  processing  at  sample  rates  above  10  MHz,  acoustic  surface  waves  are  at¬ 
tractive.  Figure  10  shows  a  possible  surface  wave  electrode  configuration.  Aluminum  elec¬ 
trodes  are  usually  deposited  on  a  piezoelectric  material  such  as  ST  quartz,  and  they  act  to 
launch  and  receive  the  acoustic  waves.  Each  transducer  is  composed  of  a  pair  of  interdigi- 
tated  electrodes  and  can  act  as  a  tap  the.  weight  of  which  is  proportional  to  the  length  of  the 
interdigitation  which  intercepts  the  acoustic  wave.  Since  the  launch  transducers  are  usually 
many  w&.elengths  long,  the  resulting  waves  are  well  collimated  and  several  filters  may  be 
implemented  simultaneously  on  a  single  substrate.  Since  both  the  registration  of  the  multi¬ 
ple  filters  and  the  weight  of  the  filter  taps  are  photo-lithographically  determined,  good 
umrormity  and  repeatability  are  possible. 
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When  it  is  desirable  to  implement  more  than  one  transform  with  the  same  hardware, 
then  programmable  filters  are  required.  These  filters  must  have  the  weighting  function 
specified  ip  advance  after  which  they  operate  as  ordinary  transversal  filters  with  one  excep¬ 
tion  -  at  the  present  time  only  binary  tap  weighting  is  possible  so  that  several  programma¬ 
ble  filters  must  be  used  simultaneously  with  a  binary  decomposition  of  the  weighting 
function. 

Surface  wave  filters  with  up  to  1  ?8  digitally  switched  taps  have  been  built.  [4,  5] . 

A  representative  programmable  ‘'ilter  is  shown  in  Figure  1  1.  It  consists  of  a  128  tap  trans¬ 
versal  filter  with  binary  taps  that  are  switched  under  the  control  of  a  shift  register.  Each  tap 
may  be  connected  to  either  of  two  busses.  The  output  signal  is  the  difference  of  the  two 
summed  bus  signals  and  is  provided  by  a  differential  transformer.  The  signal  may  be  applied 
at  either  of  transducers  A  or  B,  and  the  weight  vector  enters  at  the  terminal  labeled  “Code 
in.” 

Charge  Transfer  Devices 

For  signal  processing  at  sample  rates  up  to  10  MHz,  when  short  duration  interrup¬ 
tion  of  the  signal  processing  may  occur,  charge  transfer  devices  such  as  charge  coupled  de¬ 
vices  (CCD)  and  bucket  brigade  devices  offer  many  advantages.  These  devices  are  fabricated 
by  MOS  technology  from  silicon  and  are  small  in  size.  There  are,  however,  two  limitations, 
there  is  charge  transfer  inefficiency  and  thermally  generated  charge  combines  with  the  signal 
so  that  the  storage  time  of  the  signal  is  less  than  a  second  at  room  temperature.  Since  the 
carriers  representing  the  signal  are  moved  through  the  device  by  clock  voltages  applied  to  an 
overlying  set  of  electrodes,  the  propagation  of  the  data  signal  may  be  interrupted  by  stop¬ 
ping  the  clock.  Figure  12  adapted  from  [6]  shows  a  three  phase  CCD  transversal  filter  and 
the  method  of  sensing  the  charge.  Figure  13  shows  a  21  sample  chirp  filter  [7] . 

A  programmable  CCD  filter  is  currently  under  development  by  General  Electric  for 
the  Naval  Undersea  Center.  This  filter  will  be  built  using  a  module  size  of  64  samples  and 
will  have  a  4  MHz  sample  rate.  Its  construction  is  such  that  either  the  binary  or  the  analog 
signal  may  be  considered  as  the  input  and  the  other  signal  considered  as  the  weight.  Unlike 
the  fixed  reference  CCD  filter  previously  discussed,  this  filter  has  provision  for  cascading  of 
modules  without  degradation  of  the  signal  by  charge  transfer  inefficiency. 

Binary  Shift  Register  Crossconvolvers 

When  signal  processing  must  be  interrupted  indefinitely,  semiconductor  shift  register 
equivalents  of  programmable  transversal  filters  are  most  suitable.  Such  a  convolver  module 
(Fig.  14),  lias  been  developed  by  TRW  for  NavElex  under  the  name  of  LSI  Correlator  [8] . 
This  device  may  be  considered  as  a  bim.ry  crossconvolver  and  its  use  is  similar  to  the 


CLOCK 


(b) 

Adapted  from  IEEE  JOURNAL  OF  SOLID  STATE  CIRCUITS  (April  1973). 
in  press. 

Fig.  12.  Charge-Cou  )led  Transversal  Filter 

programmable  transversal  filters  previously  described.  However,  it  is  restricted  to  binary 
signals  and  binary  weights.  The  use  of  multiple  binary  correlation  modules  to  attain  multi¬ 
ple  level  quantization  is  described  in  [31.  Figure  15  shows  a  configuration  suitable  for 
multilevel  operation  of  one  of  the  channels.  The  LSI  correlator  is  available  in  either 
a  module  size  of  64  bits  at  a  20  MHz  clock  rate  or  a  module  size  of  8  bits  at  a  100  MHz 
clock  rate. 
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SIGNAL  "S" 


Fig.  15.  Correlator  Configuration  for  3-Bit  Signal  vs  Binary  Reference 

LINEAR  TRANSFORM  PROGRAMMABLE  SIGNAL  PROCESSING  ARITHMETIC  UNIT 

The  linear  and  bilinear  operations  most  needed  for  signal  processing  are:  matched 
filtering  or  replica  correlation,  the  discrete  Fourier  transform,  and  crosscorrelation.  These 
transforms  represent  an  excessive  computational  load  for  a  general  purpose  digital  computer 
and  a  heavy  computational  load  even  for  a  digital  computer  structured  for  signal  processing. 
The  decomposable  transform  architecture  may  be  used  to  provide  a  single  module  which  can 
perform  the  above  operations  in  a  time  proportional  to  the  data  block  length.  Such  mod¬ 
ules  may  be  configured  under  computer  control  to  permit  not  only  changes  in  task  assign¬ 
ment  but  also  an  exchange  among  data  block  length,  accuracy,  and  speed.  Figure  1 6  shows 
such  a  module  and  the  contents  of  the  local  stores  for  matched  filtering,  crosscorrelation, 
and  the  discrete  Fourier  transform. 

Serial  access  delay  line  memories  may  be  used  in  the  modular  building  blocks  for 
high  speed  non-digital  implementations.  If  interruptions  of  the  processing  are  infrequent, 
then  ultrasonic  delay  lines  are  preferred,  since  they  offer  large  storage  capacity,  convenient 
tapping,  and  low  power  dissipation.  Where  frequent  processing  interrupts  are  required, 
charge  coupled  devices  with  their  controllable  clock  rate  may  be  used.  The  corresponding 
modular  building  block  for  digital  implementations  is  a  contiguously  tapped  binary  shift 
register  crossconvolver. 
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FOR  LINEAR  TRANSFORMS 
a.  b,  c  ARE  LOCAL  STORES 


FOR  MATCH  FILTERS  AND  CROSS  CORRELATION 

a  =  1  b  =  {bj}  c  =  1 

0  <  j  <  N-1 

FOR  THE  DFT  VIA  THE  CHIRP-Z  ALGORITHM 
a  =  e-ijrn2/N  b  =  ei7rn2/N  c e-i7rn2/N 

0 < n  <  N-1  0  < n  < 2  N-2  0<i  <  N-1 

★  denotes  convolution 

Fig.  16.  Linear  Transform  Programmable  SPAU  Module 
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